source: git/Singular/walk.cc @ 35aef3

spielwiese
Last change on this file since 35aef3 was 35aef3, checked in by Stephan Oberfranz <oberfran@…>, 8 years ago
randomization fixed
  • Property mode set to 100755
File size: 231.5 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: Implementation of the Groebner walk
7*/
8
9// define if the Buchberger alg should be used
10//   to compute a reduced GB of a omega-homogenoues ideal
11// default: we use the hilbert driven algorithm.
12#define BUCHBERGER_ALG  //we use the improved Buchberger alg.
13
14//#define UPPER_BOUND //for the original "Tran" algorithm
15//#define REPRESENTATION_OF_SIGMA //if one perturbs sigma in Tran
16
17//#define TEST_OVERFLOW
18//#define CHECK_IDEAL
19//#define CHECK_IDEAL_MWALK
20
21//#define NEXT_VECTORS_CC
22//#define PRINT_VECTORS //to print vectors (sigma, tau, omega)
23
24#define INVEPS_SMALL_IN_FRACTAL  //to choose the small invers of epsilon
25#define INVEPS_SMALL_IN_MPERTVECTOR  //to choose the small invers of epsilon
26#define INVEPS_SMALL_IN_TRAN  //to choose the small invers of epsilon
27
28#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
30//                          tau doesn't stay in the correct cone
31
32//#define TIME_TEST // print the used time of each subroutine
33//#define ENDWALKS //print the size of the last omega-homogenoues Groebner basis
34
35/* includes */
36
37#include <kernel/mod2.h>
38#include <misc/intvec.h>
39#include <Singular/cntrlc.h>
40#include <misc/options.h>
41#include <omalloc/omalloc.h>
42#include <Singular/ipshell.h>
43#include <Singular/ipconv.h>
44#include <coeffs/ffields.h>
45#include <coeffs/coeffs.h>
46#include <Singular/subexpr.h>
47#include <polys/templates/p_Procs.h>
48
49#include <polys/monomials/maps.h>
50
51/* include Hilbert-function */
52#include <kernel/combinatorics/stairc.h>
53
54/** kstd2.cc */
55#include <kernel/GBEngine/kutil.h>
56#include <kernel/GBEngine/khstd.h>
57
58#include <Singular/walk.h>
59#include <kernel/polys.h>
60#include <kernel/ideals.h>
61#include <Singular/ipid.h>
62#include <Singular/tok.h>
63#include <coeffs/numbers.h>
64#include <Singular/ipid.h>
65#include <polys/monomials/ring.h>
66#include <kernel/GBEngine/kstd1.h>
67#include <polys/matpol.h>
68#include <polys/weight.h>
69#include <misc/intvec.h>
70#include <kernel/GBEngine/syz.h>
71#include <Singular/lists.h>
72#include <polys/prCopy.h>
73#include <polys/monomials/ring.h>
74//#include <polys/ext_fields/longalg.h>
75#include <polys/clapsing.h>
76
77#include <coeffs/mpr_complex.h>
78
79#include <stdio.h>
80// === Zeit & System (Holger Croeni ===
81#include <time.h>
82#include <sys/time.h>
83#include <math.h>
84#include <sys/stat.h>
85#include <unistd.h>
86#include <float.h>
87#include <misc/mylimits.h>
88#include <sys/types.h>
89
90int nstep;
91
92extern BOOLEAN ErrorCheck();
93
94extern BOOLEAN pSetm_error;
95
96void Set_Error( BOOLEAN f) { pSetm_error=f; }
97
98BOOLEAN Overflow_Error =  FALSE;
99
100clock_t xtif, xtstd, xtlift, xtred, xtnw;
101clock_t xftostd, xtextra, xftinput, to;
102
103/****************************
104 * utilities for TSet, LSet *
105 ****************************/
106inline static intset initec (int maxnr)
107{
108  return (intset)omAlloc(maxnr*sizeof(int));
109}
110
111inline static unsigned long* initsevS (int maxnr)
112{
113  return (unsigned long*)omAlloc0(maxnr*sizeof(unsigned long));
114}
115inline static int* initS_2_R (int maxnr)
116{
117  return (int*)omAlloc0(maxnr*sizeof(int));
118}
119
120/************************************
121 * construct the set s from F u {P} *
122 ************************************/
123// unused
124#if 0
125static void initSSpecialCC (ideal F, ideal Q, ideal P,kStrategy strat)
126{
127  int   i,pos;
128
129  if (Q!=NULL) i=((IDELEMS(Q)+(setmaxTinc-1))/setmaxTinc)*setmaxTinc;
130  else i=setmaxT;
131
132  strat->ecartS=initec(i);
133  strat->sevS=initsevS(i);
134  strat->S_2_R=initS_2_R(i);
135  strat->fromQ=NULL;
136  strat->Shdl=idInit(i,F->rank);
137  strat->S=strat->Shdl->m;
138
139  // - put polys into S -
140  if (Q!=NULL)
141  {
142    strat->fromQ=initec(i);
143    memset(strat->fromQ,0,i*sizeof(int));
144    for (i=0; i<IDELEMS(Q); i++)
145    {
146      if (Q->m[i]!=NULL)
147      {
148        LObject h;
149        h.p = pCopy(Q->m[i]);
150        //if (TEST_OPT_INTSTRATEGY)
151        //{
152        //  //pContent(h.p);
153        //  h.pCleardenom(); // also does a pContent
154        //}
155        //else
156        //{
157        //  h.pNorm();
158        //}
159        strat->initEcart(&h);
160        if (rHasLocalOrMixedOrdering_currRing())
161        {
162          deleteHC(&h,strat);
163        }
164        if (h.p!=NULL)
165        {
166          if (strat->sl==-1)
167            pos =0;
168          else
169          {
170            pos = posInS(strat,strat->sl,h.p,h.ecart);
171          }
172          h.sev = pGetShortExpVector(h.p);
173          h.SetpFDeg();
174          strat->enterS(h,pos,strat, strat->tl+1);
175          enterT(h, strat);
176          strat->fromQ[pos]=1;
177        }
178      }
179    }
180  }
181  //- put polys into S -
182  for (i=0; i<IDELEMS(F); i++)
183  {
184    if (F->m[i]!=NULL)
185    {
186      LObject h;
187      h.p = pCopy(F->m[i]);
188      if (rHasGlobalOrdering(currRing))
189      {
190        //h.p=redtailBba(h.p,strat->sl,strat);
191        h.p=redtailBba(h.p,strat->sl,strat);
192      }
193      else
194      {
195        deleteHC(&h,strat);
196      }
197      strat->initEcart(&h);
198      if (h.p!=NULL)
199      {
200        if (strat->sl==-1)
201          pos =0;
202        else
203          pos = posInS(strat,strat->sl,h.p,h.ecart);
204        h.sev = pGetShortExpVector(h.p);
205        strat->enterS(h,pos,strat, strat->tl+1);
206        h.length = pLength(h.p);
207        h.SetpFDeg();
208        enterT(h,strat);
209      }
210    }
211  }
212#ifdef INITSSPECIAL
213  for (i=0; i<IDELEMS(P); i++)
214  {
215    if (P->m[i]!=NULL)
216    {
217      LObject h;
218      h.p=pCopy(P->m[i]);
219      strat->initEcart(&h);
220      h.length = pLength(h.p);
221      if (TEST_OPT_INTSTRATEGY)
222      {
223        h.pCleardenom();
224      }
225      else
226      {
227        h.pNorm();
228      }
229      if(strat->sl>=0)
230      {
231        if (rHasGlobalOrdering(currRing))
232        {
233          h.p=redBba(h.p,strat->sl,strat);
234          if (h.p!=NULL)
235            h.p=redtailBba(h.p,strat->sl,strat);
236        }
237        else
238        {
239          h.p=redMora(h.p,strat->sl,strat);
240          strat->initEcart(&h);
241        }
242        if(h.p!=NULL)
243        {
244          if (TEST_OPT_INTSTRATEGY)
245          {
246            h.pCleardenom();
247          }
248          else
249          {
250            h.is_normalized = 0;
251            h.pNorm();
252          }
253          h.sev = pGetShortExpVector(h.p);
254          h.SetpFDeg();
255          pos = posInS(strat->S,strat->sl,h.p,h.ecart);
256          enterpairsSpecial(h.p,strat->sl,h.ecart,pos,strat,strat->tl+1);
257          strat->enterS(h,pos,strat, strat->tl+1);
258          enterT(h,strat);
259        }
260      }
261      else
262      {
263        h.sev = pGetShortExpVector(h.p);
264        h.SetpFDeg();
265        strat->enterS(h,0,strat, strat->tl+1);
266        enterT(h,strat);
267      }
268    }
269  }
270#endif
271}
272#endif
273
274/*****************
275 *interreduce F  *
276 *****************/
277static ideal kInterRedCC(ideal F, ideal Q)
278{
279  int j;
280  kStrategy strat = new skStrategy;
281
282//  if (TEST_OPT_PROT)
283//  {
284//    writeTime("start InterRed:");
285//    mflush();
286//  }
287  //strat->syzComp     = 0;
288  strat->kHEdgeFound = (currRing->ppNoether) != NULL;
289  strat->kNoether=pCopy((currRing->ppNoether));
290  strat->ak = id_RankFreeModule(F, currRing);
291  initBuchMoraCrit(strat);
292  strat->NotUsedAxis = (BOOLEAN *)omAlloc((currRing->N+1)*sizeof(BOOLEAN));
293  for(j=currRing->N; j>0; j--)
294  {
295    strat->NotUsedAxis[j] = TRUE;
296  }
297  strat->enterS      = enterSBba;
298  strat->posInT      = posInT0;
299  strat->initEcart   = initEcartNormal;
300  strat->sl   = -1;
301  strat->tl          = -1;
302  strat->tmax        = setmaxT;
303  strat->T           = initT();
304  strat->R           = initR();
305  strat->sevT        = initsevT();
306  if(rHasLocalOrMixedOrdering_currRing())
307  {
308    strat->honey = TRUE;
309  }
310
311  //initSCC(F,Q,strat);
312  initS(F,Q,strat);
313
314  /*
315  timetmp=clock();//22.01.02
316  initSSpecialCC(F,Q,NULL,strat);
317  tininitS=tininitS+clock()-timetmp;//22.01.02
318  */
319  if(TEST_OPT_REDSB)
320  {
321    strat->noTailReduction=FALSE;
322  }
323  updateS(TRUE,strat);
324
325  if(TEST_OPT_REDSB && TEST_OPT_INTSTRATEGY)
326  {
327    completeReduce(strat);
328  }
329  pDelete(&strat->kHEdge);
330  omFreeSize((ADDRESS)strat->T,strat->tmax*sizeof(TObject));
331  omFreeSize((ADDRESS)strat->ecartS,IDELEMS(strat->Shdl)*sizeof(int));
332  omFreeSize((ADDRESS)strat->sevS,IDELEMS(strat->Shdl)*sizeof(unsigned long));
333  omFreeSize((ADDRESS)strat->NotUsedAxis,(currRing->N+1)*sizeof(BOOLEAN));
334  omfree(strat->sevT);
335  omfree(strat->S_2_R);
336  omfree(strat->R);
337
338  if(strat->fromQ)
339  {
340    for(j=0; j<IDELEMS(strat->Shdl); j++)
341    {
342      if(strat->fromQ[j])
343      {
344        pDelete(&strat->Shdl->m[j]);
345      }
346    }
347    omFreeSize((ADDRESS)strat->fromQ,IDELEMS(strat->Shdl)*sizeof(int));
348    strat->fromQ = NULL;
349  }
350//  if (TEST_OPT_PROT)
351//  {
352//    writeTime("end Interred:");
353//    mflush();
354//  }
355  ideal shdl=strat->Shdl;
356  idSkipZeroes(shdl);
357  delete(strat);
358
359  return shdl;
360}
361
362//unused
363#if 0
364static void TimeString(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd,
365                       clock_t tlf,clock_t tred, clock_t tnw, int step)
366{
367  double totm = ((double) (clock() - tinput))/1000000;
368  double ostd,mostd, mif, mstd, mlf, mred, mnw, mxif,mxstd,mxlf,mxred,mxnw,tot;
369  // double mextra
370  Print("\n// total time = %.2f sec", totm);
371  Print("\n// tostd = %.2f sec = %.2f", ostd=((double) tostd)/1000000,
372        mostd=((((double) tostd)/1000000)/totm)*100);
373  Print("\n// tif   = %.2f sec = %.2f", ((double) tif)/1000000,
374        mif=((((double) tif)/1000000)/totm)*100);
375  Print("\n// std   = %.2f sec = %.2f", ((double) tstd)/1000000,
376        mstd=((((double) tstd)/1000000)/totm)*100);
377  Print("\n// lift  = %.2f sec = %.2f", ((double) tlf)/1000000,
378        mlf=((((double) tlf)/1000000)/totm)*100);
379  Print("\n// ired  = %.2f sec = %.2f", ((double) tred)/1000000,
380        mred=((((double) tred)/1000000)/totm)*100);
381  Print("\n// nextw = %.2f sec = %.2f", ((double) tnw)/1000000,
382        mnw=((((double) tnw)/1000000)/totm)*100);
383  PrintS("\n Time for the last step:");
384  Print("\n// xinfo = %.2f sec = %.2f", ((double) xtif)/1000000,
385        mxif=((((double) xtif)/1000000)/totm)*100);
386  Print("\n// xstd  = %.2f sec = %.2f", ((double) xtstd)/1000000,
387        mxstd=((((double) xtstd)/1000000)/totm)*100);
388  Print("\n// xlift = %.2f sec = %.2f", ((double) xtlift)/1000000,
389        mxlf=((((double) xtlift)/1000000)/totm)*100);
390  Print("\n// xired = %.2f sec = %.2f", ((double) xtred)/1000000,
391        mxred=((((double) xtred)/1000000)/totm)*100);
392  Print("\n// xnextw= %.2f sec = %.2f", ((double) xtnw)/1000000,
393        mxnw=((((double) xtnw)/1000000)/totm)*100);
394
395  tot=mostd+mif+mstd+mlf+mred+mnw+mxif+mxstd+mxlf+mxred+mxnw;
396  double res = (double) 100 - tot;
397  Print("\n// &%d&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f&%.2f(%.2f)\\ \\",
398        step, ostd, totm, mostd,mif,mstd,mlf,mred,mnw,mxif,mxstd,mxlf,mxred,mxnw,tot,res,
399        ((((double) xtextra)/1000000)/totm)*100);
400}
401#endif
402
403//unused
404#if 0
405static void TimeStringFractal(clock_t tinput, clock_t tostd, clock_t tif,clock_t tstd,
406                       clock_t textra, clock_t tlf,clock_t tred, clock_t tnw)
407{
408
409  double totm = ((double) (clock() - tinput))/1000000;
410  double ostd, mostd, mif, mstd, mextra, mlf, mred, mnw, tot, res;
411  Print("\n// total time = %.2f sec", totm);
412  Print("\n// tostd = %.2f sec = %.2f", ostd=((double) tostd)/1000000,
413        mostd=((((double) tostd)/1000000)/totm)*100);
414  Print("\n// tif   = %.2f sec = %.2f", ((double) tif)/1000000,
415        mif=((((double) tif)/1000000)/totm)*100);
416  Print("\n// std   = %.2f sec = %.2f", ((double) tstd)/1000000,
417        mstd=((((double) tstd)/1000000)/totm)*100);
418  Print("\n// xstd  = %.2f sec = %.2f", ((double) textra)/1000000,
419        mextra=((((double) textra)/1000000)/totm)*100);
420  Print("\n// lift  = %.2f sec = %.2f", ((double) tlf)/1000000,
421        mlf=((((double) tlf)/1000000)/totm)*100);
422  Print("\n// ired  = %.2f sec = %.2f", ((double) tred)/1000000,
423        mred=((((double) tred)/1000000)/totm)*100);
424  Print("\n// nextw = %.2f sec = %.2f", ((double) tnw)/1000000,
425        mnw=((((double) tnw)/1000000)/totm)*100);
426  tot = mostd+mif+mstd+mextra+mlf+mred+mnw;
427  res = (double) 100.00-tot;
428  Print("\n// &%.2f &%.2f&%.2f &%.2f &%.2f &%.2f &%.2f &%.2f &%.2f&%.2f&%.2f\\ \\ ",
429        ostd,totm,mostd,mif,mstd,mextra,mlf,mred,mnw,tot,res);
430}
431#endif
432
433//#ifdef CHECK_IDEAL_MWALK
434static void idString(ideal L, const char* st)
435{
436  int i, nL = IDELEMS(L);
437
438  Print("\n//  ideal %s =  ", st);
439  for(i=0; i<nL-1; i++)
440  {
441    Print(" %s, ", pString(L->m[i]));
442  }
443  Print(" %s;", pString(L->m[nL-1]));
444}
445//#endif
446
447#if defined(CHECK_IDEAL_MWALK) || defined(ENDWALKS)
448static void headidString(ideal L, char* st)
449{
450  int i, nL = IDELEMS(L);
451
452  Print("\n//  ideal %s =  ", st);
453  for(i=0; i<nL-1; i++)
454  {
455    Print(" %s, ", pString(pHead(L->m[i])));
456  }
457  Print(" %s;", pString(pHead(L->m[nL-1])));
458}
459#endif
460
461#if defined(CHECK_IDEAL_MWALK) || defined(ENDWALKS)
462static void idElements(ideal L, char* st)
463{
464  int i, nL = IDELEMS(L);
465  int *K=(int *)omAlloc(nL*sizeof(int));
466
467  Print("\n//  #monoms of %s =  ", st);
468  for(i=0; i<nL; i++)
469  {
470    K[i] = pLength(L->m[i]);
471  }
472  int j, nsame;
473  // int  nk=0;
474  for(i=0; i<nL; i++)
475  {
476    if(K[i]!=0)
477    {
478      nsame = 1;
479      for(j=i+1; j<nL; j++)
480      {
481        if(K[j]==K[i])
482        {
483          nsame ++;
484          K[j]=0;
485        }
486      }
487      if(nsame == 1)
488      {
489        Print("%d, ",K[i]);
490      }
491      else
492      {
493        Print("%d[%d], ", K[i], nsame);
494      }
495    }
496  }
497  omFree(K);
498}
499#endif
500
501
502static void ivString(intvec* iv, const char* ch)
503{
504  int nV = iv->length()-1;
505  Print("\n// intvec %s =  ", ch);
506
507  for(int i=0; i<nV; i++)
508  {
509    Print("%d, ", (*iv)[i]);
510  }
511  Print("%d;", (*iv)[nV]);
512}
513
514//unused
515//#if 0
516static void MivString(intvec* iva, intvec* ivb, intvec* ivc)
517{
518  int nV = iva->length()-1;
519  int i;
520  PrintS("\n//  (");
521  for(i=0; i<nV; i++)
522  {
523    Print("%d, ", (*iva)[i]);
524  }
525  Print("%d) ==> (", (*iva)[nV]);
526  for(i=0; i<nV; i++)
527  {
528    Print("%d, ", (*ivb)[i]);
529  }
530  Print("%d) := (", (*ivb)[nV]);
531
532  for(i=0; i<nV; i++)
533  {
534    Print("%d, ", (*ivc)[i]);
535  }
536  Print("%d)", (*ivc)[nV]);
537}
538//#endif
539
540/********************************************************************
541 * returns gcd of integers a and b                                  *
542 ********************************************************************/
543static inline long gcd(const long a, const long b)
544{
545  long r, p0 = a, p1 = b;
546  //assume(p0 >= 0 && p1 >= 0);
547  if(p0 < 0)
548  {
549    p0 = -p0;
550  }
551  if(p1 < 0)
552  {
553    p1 = -p1;
554  }
555  while(p1 != 0)
556  {
557    r = p0 % p1;
558    p0 = p1;
559    p1 = r;
560  }
561  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;
592}
593
594/*********************************************
595 * cancel gcd of integers zaehler and nenner *
596 *********************************************/
597static void cancel(mpz_t zaehler, mpz_t nenner)
598{
599//  assume(zaehler >= 0 && nenner > 0);
600  mpz_t g;
601  mpz_init(g);
602  mpz_gcd(g, zaehler, nenner);
603
604  mpz_div(zaehler , zaehler, g);
605  mpz_div(nenner ,  nenner, g);
606
607  mpz_clear(g);
608}
609
610//unused
611#if 0
612static int isVectorNeg(intvec* omega)
613{
614  int i;
615
616  for(i=omega->length(); i>=0; i--)
617  {
618    if((*omega)[i]<0)
619    {
620      return 1;
621    }
622  }
623  return 0;
624}
625#endif
626
627/********************************************************************
628 * compute a weight degree of a monomial p w.r.t. a weight_vector   *
629 ********************************************************************/
630static inline int MLmWeightedDegree(const poly p, intvec* weight)
631{
632  /* 2147483647 is max. integer representation in SINGULAR */
633  mpz_t sing_int;
634  mpz_init_set_ui(sing_int,  2147483647);
635
636  int i, wgrad;
637
638  mpz_t zmul;
639  mpz_init(zmul);
640  mpz_t zvec;
641  mpz_init(zvec);
642  mpz_t zsum;
643  mpz_init(zsum);
644
645  for (i=currRing->N; i>0; i--)
646  {
647    mpz_set_si(zvec, (*weight)[i-1]);
648    mpz_mul_ui(zmul, zvec, pGetExp(p, i));
649    mpz_add(zsum, zsum, zmul);
650  }
651
652  wgrad = mpz_get_ui(zsum);
653
654  if(mpz_cmp(zsum, sing_int)>0)
655  {
656    if(Overflow_Error ==  FALSE)
657    {
658      PrintLn();
659      PrintS("\n// ** OVERFLOW in \"MwalkInitialForm\": ");
660      mpz_out_str( stdout, 10, zsum);
661      PrintS(" is greater than 2147483647 (max. integer representation)");
662      Overflow_Error = TRUE;
663    }
664  }
665
666  mpz_clear(zmul);
667  mpz_clear(zvec);
668  mpz_clear(zsum);
669  mpz_clear(sing_int);
670
671  return wgrad;
672}
673
674/********************************************************************
675 * compute a weight degree of a polynomial p w.r.t. a weight_vector *
676 ********************************************************************/
677static inline int MwalkWeightDegree(poly p, intvec* weight_vector)
678{
679  assume(weight_vector->length() >= currRing->N);
680  int max = 0, maxtemp;
681
682  while(p != NULL)
683  {
684    maxtemp = MLmWeightedDegree(p, weight_vector);
685    pIter(p);
686
687    if (maxtemp > max)
688    {
689      max = maxtemp;
690    }
691  }
692  return max;
693}
694
695
696/********************************************************************
697 * compute a weight degree of a monomial p w.r.t. a weight_vector   *
698 ********************************************************************/
699static  void  MLmWeightedDegree_gmp(mpz_t result, const poly p, intvec* weight)
700{
701  /* 2147483647 is max. integer representation in SINGULAR */
702  mpz_t sing_int;
703  mpz_init_set_ui(sing_int,  2147483647);
704
705  int i;
706
707  mpz_t zmul;
708  mpz_init(zmul);
709  mpz_t zvec;
710  mpz_init(zvec);
711  mpz_t ztmp;
712  mpz_init(ztmp);
713
714  for (i=currRing->N; i>0; i--)
715  {
716    mpz_set_si(zvec, (*weight)[i-1]);
717    mpz_mul_ui(zmul, zvec, pGetExp(p, i));
718    mpz_add(ztmp, ztmp, zmul);
719  }
720  mpz_init_set(result, ztmp);
721  mpz_clear(ztmp);
722  mpz_clear(sing_int);
723  mpz_clear(zvec);
724  mpz_clear(zmul);
725}
726
727
728/*****************************************************************************
729 * return an initial form of the polynom g w.r.t. a weight vector curr_weight *
730 *****************************************************************************/
731static poly MpolyInitialForm(poly g, intvec* curr_weight)
732{
733  if(g == NULL)
734  {
735    return NULL;
736  }
737  mpz_t max; mpz_init(max);
738  mpz_t maxtmp; mpz_init(maxtmp);
739
740  poly hg, in_w_g = NULL;
741
742  while(g != NULL)
743  {
744    hg = g;
745    pIter(g);
746    MLmWeightedDegree_gmp(maxtmp, hg, curr_weight);
747
748    if(mpz_cmp(maxtmp, max)>0)
749    {
750      mpz_init_set(max, maxtmp);
751      pDelete(&in_w_g);
752      in_w_g = pHead(hg);
753    }
754    else
755    {
756      if(mpz_cmp(maxtmp, max)==0)
757      {
758        in_w_g = pAdd(in_w_g, pHead(hg));
759      }
760    }
761  }
762  return in_w_g;
763}
764
765/************************************************************************
766 * compute the initial form of an ideal <G> w.r.t. a weight vector iva  *
767 ************************************************************************/
768ideal MwalkInitialForm(ideal G, intvec* ivw)
769{
770  BOOLEAN nError =  Overflow_Error;
771  Overflow_Error = FALSE;
772
773  int i, nG = IDELEMS(G);
774  ideal Gomega = idInit(nG, 1);
775
776  for(i=nG-1; i>=0; i--)
777  {
778    Gomega->m[i] = MpolyInitialForm(G->m[i], ivw);
779  }
780  if(Overflow_Error == FALSE)
781  {
782    Overflow_Error = nError;
783  }
784  return Gomega;
785}
786
787/************************************************************************
788 * test whether the weight vector iv is in the cone of the ideal G      *
789 *     i.e. test whether in(in_w(g)) = in(g) for all g in G             *
790 ************************************************************************/
791
792static int test_w_in_ConeCC(ideal G, intvec* iv)
793{
794  if(G->m[0] == NULL)
795  {
796    PrintS("//** the result may be WRONG, i.e. 0!!\n");
797    return 0;
798  }
799
800  BOOLEAN nError =  Overflow_Error;
801  Overflow_Error = FALSE;
802
803  int i, nG = IDELEMS(G);
804  poly mi, gi;
805
806  for(i=nG-1; i>=0; i--)
807  {
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));
812    if(mi == NULL)
813    {
814      pDelete(&mi);
815      if(Overflow_Error == FALSE)
816      {
817        Overflow_Error = nError;
818      }
819      return 0;
820    }
821    if(!pLmEqual(mi, gi))
822    {
823      pDelete(&mi);
824      if(Overflow_Error == FALSE)
825      {
826        Overflow_Error = nError;
827      }
828      return 0;
829    }
830    pDelete(&mi);
831  }
832
833  if(Overflow_Error == FALSE)
834  {
835    Overflow_Error = nError;
836  }
837  return 1;
838}
839
840/***************************************************
841 * compute a least common multiple of two integers *
842 ***************************************************/
843static inline long Mlcm(long &i1, long &i2)
844{
845  long temp = gcd(i1, i2);
846  return ((i1 / temp)* i2);
847}
848
849
850/***************************************************
851 * return  the dot product of two intvecs a and b  *
852 ***************************************************/
853static inline long  MivDotProduct(intvec* a, intvec* b)
854{
855  assume( a->length() ==  b->length());
856  int i, n = a->length();
857  long result = 0;
858
859  for(i=n-1; i>=0; i--)
860    {
861    result += (*a)[i] * (*b)[i];
862    }
863  return result;
864}
865
866/*****************************************************
867 * Substract two given intvecs componentwise         *
868 *****************************************************/
869static intvec* MivSub(intvec* a, intvec* b)
870{
871  assume( a->length() ==  b->length());
872  int i, n = a->length();
873  intvec* result = new intvec(n);
874
875  for(i=n-1; i>=0; i--)
876  {
877    (*result)[i] = (*a)[i] - (*b)[i];
878  }
879  return result;
880}
881
882/*****************************************************
883 * return the "intvec" lead exponent of a polynomial *
884 *****************************************************/
885static intvec* MExpPol(poly f)
886{
887  int i, nR = currRing->N;
888  intvec* result = new intvec(nR);
889
890  for(i=nR-1; i>=0; i--)
891  {
892    (*result)[i] = pGetExp(f,i+1);
893  }
894  return result;
895}
896
897/*****************************************************
898 * Compare two given intvecs and return 1, if they   *
899 * are the same, otherwise 0                         *
900 *****************************************************/
901int MivSame(intvec* u , intvec* v)
902{
903  assume(u->length() == v->length());
904
905  int i, niv = u->length();
906
907  for (i=0; i<niv; i++)
908  {
909    if ((*u)[i] != (*v)[i])
910    {
911      return 0;
912    }
913  }
914  return 1;
915}
916
917/******************************************************
918 * Compare 3 given intvecs and return 0, if the first *
919 * and the second are the same. Return 1, if the      *
920 * the second and the third are the same, otherwise 2 *
921 ******************************************************/
922int M3ivSame(intvec* temp, intvec* u , intvec* v)
923{
924  assume(temp->length() == u->length() && u->length() == v->length());
925
926  if((MivSame(temp, u)) == 1)
927  {
928    return 0;
929  }
930  if((MivSame(temp, v)) == 1)
931  {
932    return 1;
933  }
934  return 2;
935}
936
937/*****************************************************
938 * compute a Groebner basis of an ideal              *
939 *****************************************************/
940static ideal MstdCC(ideal G)
941{
942  BITSET save1,save2;
943  SI_SAVE_OPT(save1,save2);
944  si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
945  ideal G1 = kStd(G, NULL, testHomog, NULL);
946  SI_RESTORE_OPT(save1,save2);
947
948  idSkipZeroes(G1);
949  return G1;
950}
951
952/*****************************************************
953 * compute a Groebner basis of an homogeneous ideal  *
954 *****************************************************/
955static ideal MstdhomCC(ideal G)
956{
957  BITSET save1,save2;
958  SI_SAVE_OPT(save1,save2);
959  si_opt_1|=(Sy_bit(OPT_REDTAIL)|Sy_bit(OPT_REDSB));
960  ideal G1 = kStd(G, NULL, isHomog, NULL);
961  SI_RESTORE_OPT(save1,save2);
962
963  idSkipZeroes(G1);
964  return G1;
965}
966
967
968/*****************************************************************************
969* create a weight matrix order as intvec of an extra weight vector (a(iv),lp)*
970******************************************************************************/
971intvec* MivMatrixOrder(intvec* iv)
972{
973  int i, nR = iv->length();
974
975  intvec* ivm = new intvec(nR*nR);
976
977  for(i=0; i<nR; i++)
978  {
979    (*ivm)[i] = (*iv)[i];
980  }
981  for(i=1; i<nR; i++)
982  {
983    (*ivm)[i*nR+i-1] = 1;
984  }
985  return ivm;
986}
987
988/*********************************************************************************
989* create a weight matrix order as intvec of an extra weight vector (a(iv),M(iw)) *
990**********************************************************************************/
991intvec* MivMatrixOrderRefine(intvec* iv, intvec* iw)
992{
993  assume((iv->length())*(iv->length()) == iw->length());
994  int i,j, nR = iv->length();
995 
996  intvec* ivm = new intvec(nR*nR);
997
998  for(i=0; i<nR; i++)
999  {
1000    (*ivm)[i] = (*iv)[i];
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    }
1008  }
1009  return ivm;
1010}
1011
1012/*******************************
1013 * return intvec = (1, ..., 1) *
1014 *******************************/
1015intvec* Mivdp(int nR)
1016{
1017  int i;
1018  intvec* ivm = new intvec(nR);
1019
1020  for(i=nR-1; i>=0; i--)
1021  {
1022    (*ivm)[i] = 1;
1023  }
1024  return ivm;
1025}
1026
1027/**********************************
1028 * return intvvec = (1,0, ..., 0) *
1029 **********************************/
1030intvec* Mivlp(int nR)
1031{
1032  intvec* ivm = new intvec(nR);
1033  (*ivm)[0] = 1;
1034
1035  return ivm;
1036}
1037
1038//unused
1039/*****************************************************************************
1040 * print the max total degree and the max coefficient of G                   *
1041 *****************************************************************************/
1042#if 0
1043static void checkComplexity(ideal G, char* cG)
1044{
1045  int nV = currRing->N;
1046  int nG = IDELEMS(G);
1047  intvec* ivUnit = Mivdp(nV);
1048  int i, tmpdeg, maxdeg=0;
1049  number tmpcoeff , maxcoeff=currRing->cf->nNULL;
1050  poly p;
1051  for(i=nG-1; i>=0; i--)
1052  {
1053    tmpdeg = MwalkWeightDegree(G->m[i], ivUnit);
1054    if(tmpdeg > maxdeg )
1055    {
1056      maxdeg = tmpdeg;
1057    }
1058  }
1059
1060  for(i=nG-1; i>=0; i--)
1061  {
1062    p = pCopy(G->m[i]);
1063    while(p != NULL)
1064    {
1065      //tmpcoeff = pGetCoeff(pHead(p));
1066      tmpcoeff = pGetCoeff(p);
1067      if(nGreater(tmpcoeff,maxcoeff))
1068      {
1069         maxcoeff = nCopy(tmpcoeff);
1070      }
1071      pIter(p);
1072    }
1073    pDelete(&p);
1074  }
1075  p = pNSet(maxcoeff);
1076  char* pStr = pString(p);
1077  delete ivUnit;
1078  Print("// max total degree of %s = %d\n",cG, maxdeg);
1079  Print("// max coefficient of %s  = %s", cG, pStr);//ing(p));
1080  Print(" which consists of %d digits", (int)strlen(pStr));
1081  PrintLn();
1082}
1083#endif
1084
1085/*****************************************************************************
1086* If target_ord = intmat(A1, ..., An) then calculate the perturbation        *
1087* vectors                                                                    *
1088*   tau_p_dep = inveps^(p_deg-1)*A1 + inveps^(p_deg-2)*A2 +... + A_p_deg     *
1089* where                                                                      *
1090*      inveps > totaldegree(G)*(max(A2)+...+max(A_p_deg))                    *
1091* intmat target_ord is an integer order matrix of the monomial ordering of   *
1092* basering.                                                                  *
1093* This programm computes a perturbated vector with a p_deg perturbation      *
1094* degree which smaller than the numbers of variables                         *
1095******************************************************************************/
1096intvec* MPertVectors(ideal G, intvec* ivtarget, int pdeg)
1097{
1098  // ivtarget is a matrix order of a degree reverse lex. order
1099  int nV = currRing->N;
1100  //assume(pdeg <= nV && pdeg >= 0);
1101
1102  int i, j, nG = IDELEMS(G);
1103  intvec* v_null =  new intvec(nV);
1104
1105
1106  // Check that the perturbed degree is valid
1107  if(pdeg > nV || pdeg <= 0)
1108  {
1109    WerrorS("//** The perturbed degree is wrong!!");
1110    return v_null;
1111  }
1112  delete v_null;
1113
1114  if(pdeg == 1)
1115  {
1116    return ivtarget;
1117  }
1118  mpz_t *pert_vector = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
1119  //mpz_t *pert_vector1 = (mpz_t*)omAlloc(nV*sizeof(mpz_t));
1120
1121  for(i=0; i<nV; i++)
1122  {
1123    mpz_init_set_si(pert_vector[i], (*ivtarget)[i]);
1124   // mpz_init_set_si(pert_vector1[i], (*ivtarget)[i]);
1125  }
1126  // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
1127  // where the Ai are the i-te rows of the matrix target_ord.
1128  int ntemp, maxAi, maxA=0;
1129  for(i=1; i<pdeg; i++)
1130  {
1131    maxAi = (*ivtarget)[i*nV];
1132    if(maxAi<0)
1133    {
1134      maxAi = -maxAi;
1135    }
1136    for(j=i*nV+1; j<(i+1)*nV; j++)
1137    {
1138      ntemp = (*ivtarget)[j];
1139      if(ntemp < 0)
1140      {
1141        ntemp = -ntemp;
1142      }
1143      if(ntemp > maxAi)
1144      {
1145        maxAi = ntemp;
1146      }
1147    }
1148    maxA += maxAi;
1149  }
1150
1151  // Calculate inveps = 1/eps, where 1/eps > totaldeg(p)*max1 for all p in G.
1152
1153  intvec* ivUnit = Mivdp(nV);
1154
1155  mpz_t tot_deg; mpz_init(tot_deg);
1156  mpz_t maxdeg; mpz_init(maxdeg);
1157  mpz_t inveps; mpz_init(inveps);
1158
1159
1160  for(i=nG-1; i>=0; i--)
1161  {
1162     mpz_set_ui(maxdeg, MwalkWeightDegree(G->m[i], ivUnit));
1163     if (mpz_cmp(maxdeg,  tot_deg) > 0 )
1164     {
1165       mpz_set(tot_deg, maxdeg);
1166     }
1167  }
1168
1169  delete ivUnit;
1170  mpz_mul_ui(inveps, tot_deg, maxA);
1171  mpz_add_ui(inveps, inveps, 1);
1172
1173
1174  // takes  "small" inveps
1175#ifdef INVEPS_SMALL_IN_MPERTVECTOR
1176  if(mpz_cmp_ui(inveps, pdeg)>0 && pdeg > 3)
1177  {
1178    //  Print("\n// choose the\"small\" inverse epsilon := %d / %d = ", mpz_get_si(inveps), pdeg);
1179    mpz_fdiv_q_ui(inveps, inveps, pdeg);
1180    // mpz_out_str(stdout, 10, inveps);
1181  }
1182#else
1183  // PrintS("\n// the \"big\" inverse epsilon: ");
1184  mpz_out_str(stdout, 10, inveps);
1185#endif
1186
1187  // pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg,
1188  // pert_vector := A1
1189  for( i=1; i < pdeg; i++ )
1190  {
1191    for(j=0; j<nV; j++)
1192    {
1193      mpz_mul(pert_vector[j], pert_vector[j], inveps);
1194      if((*ivtarget)[i*nV+j]<0)
1195      {
1196        mpz_sub_ui(pert_vector[j], pert_vector[j],-(*ivtarget)[i*nV+j]);
1197      }
1198      else
1199      {
1200        mpz_add_ui(pert_vector[j], pert_vector[j],(*ivtarget)[i*nV+j]);
1201      }
1202    }
1203  }
1204  mpz_t ztemp;
1205  mpz_init(ztemp);
1206  mpz_set(ztemp, pert_vector[0]);
1207  for(i=1; i<nV; i++)
1208  {
1209    mpz_gcd(ztemp, ztemp, pert_vector[i]);
1210    if(mpz_cmp_si(ztemp, 1)  == 0)
1211    {
1212      break;
1213    }
1214  }
1215  if(mpz_cmp_si(ztemp, 1) != 0)
1216  {
1217    for(i=0; i<nV; i++)
1218    {
1219      mpz_divexact(pert_vector[i], pert_vector[i], ztemp);
1220    }
1221  }
1222
1223  intvec *pert_vector1= new intvec(nV);
1224  j = 0;
1225  for(i=0; i<nV; i++)
1226  {
1227    (* pert_vector1)[i] = mpz_get_si(pert_vector[i]);
1228    (* pert_vector1)[i] = 0.1*(* pert_vector1)[i];
1229    (* pert_vector1)[i] = floor((* pert_vector1)[i] + 0.5);
1230    if((* pert_vector1)[i] == 0)
1231    {
1232      j++;
1233    }
1234  }
1235  if(j > nV - 1)
1236  {
1237    // Print("\n//  MPertVectors: geaenderter vector gleich Null! \n");
1238    delete pert_vector1;
1239    goto CHECK_OVERFLOW;
1240  }
1241
1242// check that the perturbed weight vector lies in the Groebner cone
1243  if(test_w_in_ConeCC(G,pert_vector1) != 0)
1244  {
1245    // Print("\n//  MPertVectors: geaenderter vector liegt in Groebnerkegel! \n");
1246    for(i=0; i<nV; i++)
1247    {
1248      mpz_set_si(pert_vector[i], (*pert_vector1)[i]);
1249    }
1250  }
1251  else
1252  {
1253    //Print("\n// MpertVectors: geaenderter vector liegt nicht in Groebnerkegel! \n");
1254  }
1255  delete pert_vector1;
1256
1257  CHECK_OVERFLOW:
1258  intvec* result = new intvec(nV);
1259
1260  /* 2147483647 is max. integer representation in SINGULAR */
1261  mpz_t sing_int;
1262  mpz_init_set_ui(sing_int,  2147483647);
1263
1264  int ntrue=0;
1265  for(i=0; i<nV; i++)
1266  {
1267    (*result)[i] = mpz_get_si(pert_vector[i]);
1268    if(mpz_cmp(pert_vector[i], sing_int)>=0)
1269    {
1270      ntrue++;
1271      if(Overflow_Error == FALSE)
1272      {
1273        Overflow_Error = TRUE;
1274        PrintS("\n// ** OVERFLOW in \"MPertvectors\": ");
1275        mpz_out_str( stdout, 10, pert_vector[i]);
1276        PrintS(" is greater than 2147483647 (max. integer representation)");
1277        Print("\n//  So vector[%d] := %d is wrong!!", i+1, (*result)[i]);
1278      }
1279    }
1280  }
1281
1282  if(Overflow_Error == TRUE)
1283  {
1284    ivString(result, "pert_vector");
1285    Print("\n// %d element(s) of it is overflow!!", ntrue);
1286  }
1287
1288  mpz_clear(ztemp);
1289  mpz_clear(sing_int);
1290  omFree(pert_vector);
1291  //omFree(pert_vector1);
1292  mpz_clear(tot_deg);
1293  mpz_clear(maxdeg);
1294  mpz_clear(inveps);
1295
1296  rComplete(currRing);
1297  for(j=0; j<IDELEMS(G); j++)
1298  {
1299    poly p=G->m[j];
1300    while(p!=NULL)
1301    {
1302      p_Setm(p,currRing); pIter(p);
1303    }
1304  }
1305  return result;
1306}
1307
1308/*****************************************************************************
1309 * The following procedure returns                                           *
1310 *     Pert(A1) = 1/eps^(pdeg-1)*A_1 + 1/eps^(pdeg-2)*A_2+...+A_pdeg,        *
1311 * where the A_i are the i-th rows of the matrix target_ord and              *
1312 *     1/eps > deg(p)*(max(A_2) + max(A_3)+...+max(A_pdeg))                  *
1313 *****************************************************************************/
1314intvec* MPertVectorslp(ideal G, intvec* ivtarget, int pdeg)
1315{
1316  // ivtarget is a matrix order of the lex. order
1317  int nV = currRing->N;
1318  //assume(pdeg <= nV && pdeg >= 0);
1319
1320  int i, j, nG = IDELEMS(G);
1321  intvec* pert_vector =  new intvec(nV);
1322
1323  //Checking that the perturbated degree is valid
1324  if(pdeg > nV || pdeg <= 0)
1325  {
1326    WerrorS("//** The perturbed degree is wrong!!");
1327    return pert_vector;
1328  }
1329  for(i=0; i<nV; i++)
1330  {
1331    (*pert_vector)[i]=(*ivtarget)[i];
1332  }
1333  if(pdeg == 1)
1334  {
1335    return pert_vector;
1336  }
1337  // Calculate max1 = Max(A2)+Max(A3)+...+Max(Apdeg),
1338  // where the Ai are the i-te rows of the matrix target_ord.
1339  int ntemp, maxAi, maxA=0;
1340  for(i=1; i<pdeg; i++)
1341  {
1342    maxAi = (*ivtarget)[i*nV];
1343    for(j=i*nV+1; j<(i+1)*nV; j++)
1344    {
1345      ntemp = (*ivtarget)[j];
1346      if(ntemp > maxAi)
1347      {
1348        maxAi = ntemp;
1349      }
1350    }
1351    maxA += maxAi;
1352  }
1353
1354  // Calculate inveps := 1/eps, where 1/eps > deg(p)*max1 for all p in G.
1355  int inveps, tot_deg = 0, maxdeg;
1356
1357  intvec* ivUnit = Mivdp(nV);//19.02
1358  for(i=nG-1; i>=0; i--)
1359  {
1360    // maxdeg = pTotaldegree(G->m[i], currRing); //it's wrong for ex1,2,rose
1361    maxdeg = MwalkWeightDegree(G->m[i], ivUnit);
1362    if (maxdeg > tot_deg )
1363    {
1364      tot_deg = maxdeg;
1365    }
1366  }
1367  delete ivUnit;
1368
1369  inveps = (tot_deg * maxA) + 1;
1370
1371#ifdef INVEPS_SMALL_IN_FRACTAL
1372  //  Print("\n// choose the\"small\" inverse epsilon := %d / %d = ", inveps, pdeg);
1373  if(inveps > pdeg && pdeg > 3)
1374  {
1375    inveps = inveps / pdeg;
1376  }
1377  // Print(" %d", inveps);
1378#else
1379  PrintS("\n// the \"big\" inverse epsilon %d", inveps);
1380#endif
1381
1382  // Pert(A1) = inveps^(pdeg-1)*A1 + inveps^(pdeg-2)*A2+...+A_pdeg
1383  for ( i=1; i < pdeg; i++ )
1384  {
1385    for(j=0; j<nV; j++)
1386    {
1387      (*pert_vector)[j] = inveps*((*pert_vector)[j]) + (*ivtarget)[i*nV+j];
1388    }
1389  }
1390
1391  int temp = (*pert_vector)[0];
1392  for(i=1; i<nV; i++)
1393  {
1394    temp = gcd(temp, (*pert_vector)[i]);
1395    if(temp == 1)
1396    {
1397      break;
1398    }
1399  }
1400  if(temp != 1)
1401  {
1402    for(i=0; i<nV; i++)
1403    {
1404      (*pert_vector)[i] = (*pert_vector)[i] / temp;
1405    }
1406  }
1407
1408  intvec* result = pert_vector;
1409  delete pert_vector;
1410  return result;
1411}
1412
1413/*****************************************************************************
1414 * define a lexicographic order matrix as intvec                             *
1415 *****************************************************************************/
1416intvec* MivMatrixOrderlp(int nV)
1417{
1418  int i;
1419  intvec* ivM = new intvec(nV*nV);
1420
1421  for(i=0; i<nV; i++)
1422  {
1423    (*ivM)[i*nV + i] = 1;
1424  }
1425  return(ivM);
1426}
1427
1428
1429/*****************************************************************************
1430 * define a reverse lexicographic order (dp) matrix as intvec                *
1431 *****************************************************************************/
1432intvec* MivMatrixOrderdp(int nV)
1433{
1434  int i;
1435  intvec* ivM = new intvec(nV*nV);
1436
1437  for(i=0; i<nV; i++)
1438  {
1439    (*ivM)[i] = 1;
1440  }
1441  for(i=1; i<nV; i++)
1442  {
1443    (*ivM)[(i+1)*nV - i] = -1;
1444  }
1445  return(ivM);
1446}
1447
1448/*****************************************************************************
1449 * creates an intvec of the monomial order Wp(ivstart)                       *
1450 *****************************************************************************/
1451intvec* MivWeightOrderlp(intvec* ivstart)
1452{
1453  int i;
1454  int nV = ivstart->length();
1455  intvec* ivM = new intvec(nV*nV);
1456
1457  for(i=0; i<nV; i++)
1458  {
1459    (*ivM)[i] = (*ivstart)[i];
1460  }
1461  for(i=1; i<nV; i++)
1462  {
1463    (*ivM)[i*nV + i-1] = 1;
1464  }
1465  return(ivM);
1466}
1467
1468/*****************************************************************************
1469 * creates an intvec of the monomial order dp(ivstart)                       *
1470 *****************************************************************************/
1471intvec* MivWeightOrderdp(intvec* ivstart)
1472{
1473  int i;
1474  int nV = ivstart->length();
1475  intvec* ivM = new intvec(nV*nV);
1476
1477  for(i=0; i<nV; i++)
1478  {
1479    (*ivM)[i] = (*ivstart)[i];
1480  }
1481  for(i=0; i<nV; i++)
1482  {
1483    (*ivM)[nV+i] = 1;
1484  }
1485  for(i=2; i<nV; i++)
1486  {
1487    (*ivM)[(i+1)*nV - i] = -1;
1488  }
1489  return(ivM);
1490}
1491
1492//unused
1493#if 0
1494static intvec* MatrixOrderdp(int nV)
1495{
1496  int i;
1497  intvec* ivM = new intvec(nV*nV);
1498
1499  for(i=0; i<nV; i++)
1500  {
1501    (*ivM)[i] = 1;
1502  }
1503  for(i=1; i<nV; i++)
1504  {
1505    (*ivM)[(i+1)*nV - i] = -1;
1506  }
1507  return(ivM);
1508}
1509#endif
1510
1511intvec* MivUnit(int nV)
1512{
1513  int i;
1514  intvec* ivM = new intvec(nV);
1515  for(i=nV-1; i>=0; i--)
1516  {
1517    (*ivM)[i] = 1;
1518  }
1519  return(ivM);
1520}
1521
1522
1523/************************************************************************
1524*  compute a perturbed weight vector of a matrix order w.r.t. an ideal  *
1525*************************************************************************/
1526int Xnlev;
1527intvec* Mfpertvector(ideal G, intvec* ivtarget)
1528{
1529  int i, j, nG = IDELEMS(G);
1530  int nV = currRing->N;
1531  int niv = nV*nV;
1532
1533
1534  // Calculate maxA = Max(A2) + Max(A3) + ... + Max(AnV),
1535  // where the Ai are the i-te rows of the matrix 'targer_ord'.
1536  int ntemp, maxAi, maxA=0;
1537  for(i=1; i<nV; i++)
1538  {
1539    maxAi = (*ivtarget)[i*nV];
1540    if(maxAi<0)
1541    {
1542      maxAi = -maxAi;
1543    }
1544    for(j=i*nV+1; j<(i+1)*nV; j++)
1545    {
1546      ntemp = (*ivtarget)[j];
1547      if(ntemp < 0)
1548      {
1549        ntemp = -ntemp;
1550      }
1551      if(ntemp > maxAi)
1552      {
1553        maxAi = ntemp;
1554      }
1555    }
1556    maxA = maxA + maxAi;
1557  }
1558  intvec* ivUnit = Mivdp(nV);
1559
1560  // Calculate inveps = 1/eps, where 1/eps > deg(p)*maxA for all p in G.
1561  mpz_t tot_deg; mpz_init(tot_deg);
1562  mpz_t maxdeg; mpz_init(maxdeg);
1563  mpz_t inveps; mpz_init(inveps);
1564
1565
1566  for(i=nG-1; i>=0; i--)
1567  {
1568    mpz_set_ui(maxdeg, MwalkWeightDegree(G->m[i], ivUnit));
1569    if (mpz_cmp(maxdeg,  tot_deg) > 0 )
1570    {
1571      mpz_set(tot_deg, maxdeg);
1572    }
1573  }
1574
1575  delete ivUnit;
1576  //inveps = (tot_deg * maxA) + 1;
1577  mpz_mul_ui(inveps, tot_deg, maxA);
1578  mpz_add_ui(inveps, inveps, 1);
1579
1580  // takes  "small" inveps
1581#ifdef INVEPS_SMALL_IN_FRACTAL
1582  if(mpz_cmp_ui(inveps, nV)>0 && nV > 3)
1583  {
1584    mpz_cdiv_q_ui(inveps, inveps, nV);
1585  }
1586  //PrintS("\n// choose the \"small\" inverse epsilon!");
1587#endif
1588
1589  // PrintLn();  mpz_out_str(stdout, 10, inveps);
1590
1591  // Calculate the perturbed target orders:
1592  mpz_t *ivtemp=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
1593  mpz_t *pert_vector=(mpz_t *)omAlloc(niv*sizeof(mpz_t));
1594
1595  for(i=0; i < nV; i++)
1596  {
1597    mpz_init_set_si(ivtemp[i], (*ivtarget)[i]);
1598    mpz_init_set_si(pert_vector[i], (*ivtarget)[i]);
1599  }
1600
1601  mpz_t ztmp; mpz_init(ztmp);
1602  // BOOLEAN isneg = FALSE;
1603
1604  for(i=1; i<nV; i++)
1605  {
1606    for(j=0; j<nV; j++)
1607    {
1608      mpz_mul(ztmp, inveps, ivtemp[j]);
1609      if((*ivtarget)[i*nV+j]<0)
1610      {
1611        mpz_sub_ui(ivtemp[j], ztmp, -(*ivtarget)[i*nV+j]);
1612      }
1613      else
1614      {
1615        mpz_add_ui(ivtemp[j], ztmp,(*ivtarget)[i*nV+j]);
1616      }
1617    }
1618
1619    for(j=0; j<nV; j++)
1620      {
1621      mpz_init_set(pert_vector[i*nV+j],ivtemp[j]);
1622      }
1623  }
1624
1625  /* 2147483647 is max. integer representation in SINGULAR */
1626  mpz_t sing_int;
1627  mpz_init_set_ui(sing_int,  2147483647);
1628
1629  intvec* result = new intvec(niv);
1630  intvec* result1 = new intvec(niv);
1631  BOOLEAN nflow = FALSE;
1632
1633  // computes gcd
1634  mpz_set(ztmp, pert_vector[0]);
1635  for(i=0; i<niv; i++)
1636  {
1637    mpz_gcd(ztmp, ztmp, pert_vector[i]);
1638    if(mpz_cmp_si(ztmp, 1)==0)
1639    {
1640      break;
1641    }
1642  }
1643
1644  for(i=0; i<niv; i++)
1645  {
1646    mpz_divexact(pert_vector[i], pert_vector[i], ztmp);
1647    (* result)[i] = mpz_get_si(pert_vector[i]);
1648  }
1649/*
1650  j = 0;
1651  for(i=0; i<niv; i++)
1652  {
1653    (* result1)[i] = mpz_get_si(pert_vector[i]);
1654    (* result1)[i] = 0.1*(* result1)[i];
1655    (* result1)[i] = floor((* result1)[i] + 0.5);
1656    if((* result1)[i] == 0)
1657    {
1658      j++;
1659    }
1660  }
1661  if(j > niv - 1)
1662  {
1663    // Print("\n//  MfPertwalk: geaenderter vector gleich Null! \n");
1664    delete result1;
1665    goto CHECK_OVERFLOW;
1666  }
1667/*
1668// check that the perturbed weight vector lies in the Groebner cone
1669  Print("\n========================================\n//** MfPertvector: test in Cone.\n");
1670  if(test_w_in_ConeCC(G,result1) != 0)
1671  {
1672    // Print("\n//  MfPertwalk: geaenderter vector liegt in Groebnerkegel! \n");
1673    delete result;
1674    result = result1;
1675    for(i=0; i<nV; i++)
1676    {
1677      mpz_set_si(pert_vector[i], (*result1)[i]);
1678    }
1679  }
1680  else
1681  {
1682    delete result1;
1683    // Print("\n// Mfpertwalk: geaenderter vector liegt nicht in Groebnerkegel! \n");
1684  }
1685  Print("\n ========================================\n");*/
1686  CHECK_OVERFLOW:
1687
1688  for(i=0; i<niv; i++)
1689  {
1690    if(mpz_cmp(pert_vector[i], sing_int)>0)
1691    {
1692      if(nflow == FALSE)
1693      {
1694        Xnlev = i / nV;
1695        nflow = TRUE;
1696        Overflow_Error = TRUE;
1697        Print("\n// Xlev = %d and the %d-th element is", Xnlev,  i+1);
1698        PrintS("\n// ** OVERFLOW in \"Mfpertvector\": ");
1699        mpz_out_str( stdout, 10, pert_vector[i]);
1700        PrintS(" is greater than 2147483647 (max. integer representation)");
1701        Print("\n//  So vector[%d] := %d is wrong!!", i+1, (*result)[i]);
1702      }
1703    }
1704  }
1705  if(Overflow_Error == TRUE)
1706  {
1707    ivString(result, "new_vector");
1708  }
1709  omFree(pert_vector);
1710  omFree(ivtemp);
1711  mpz_clear(ztmp);
1712  mpz_clear(tot_deg);
1713  mpz_clear(maxdeg);
1714  mpz_clear(inveps);
1715  mpz_clear(sing_int);
1716
1717  rComplete(currRing);
1718  for(j=0; j<IDELEMS(G); j++)
1719  {
1720    poly p=G->m[j];
1721    while(p!=NULL)
1722    {
1723      p_Setm(p,currRing); pIter(p);
1724    }
1725  }
1726  return result;
1727}
1728
1729/****************************************************************
1730 * Multiplication of two ideals element by element              *
1731 * i.e. Let be A := (a_i) and B := (b_i), return C := (a_i*b_i) *
1732 *  destroy A, keeps B                                          *
1733 ****************************************************************/
1734static ideal MidMult(ideal A, ideal B)
1735{
1736  int mA = IDELEMS(A), mB = IDELEMS(B);
1737
1738  if(A==NULL || B==NULL)
1739  {
1740    return NULL;
1741  }
1742  if(mB < mA)
1743  {
1744    mA = mB;
1745  }
1746  ideal result = idInit(mA, 1);
1747
1748  int i, k=0;
1749  for(i=0; i<mA; i++)
1750    {
1751      result->m[k] = pMult(A->m[i], pCopy(B->m[i]));
1752      A->m[i]=NULL;
1753      if (result->m[k]!=NULL)
1754      {
1755        k++;
1756      }
1757    }
1758
1759  idDelete(&A);
1760  idSkipZeroes(result);
1761  return result;
1762}
1763
1764/*********************************************************************
1765 * G is a red. Groebner basis w.r.t. <_1                             *
1766 * Gomega is an initial form ideal of <G> w.r.t. a weight vector w   *
1767 * M is a subideal of <Gomega> and M selft is a red. Groebner basis  *
1768 *    of the ideal <Gomega> w.r.t. <_w                               *
1769 * Let m_i = h1.gw1 + ... + hs.gws for each m_i in M; gwi in Gomega  *
1770 * return F with n(F) = n(M) and f_i = h1.g1 + ... + hs.gs for each i*
1771 ********************************************************************/
1772static ideal MLifttwoIdeal(ideal Gw, ideal M, ideal G)
1773{
1774  ideal Mtmp = idLift(Gw, M, NULL, FALSE, TRUE, TRUE, NULL);
1775
1776  // If Gw is a GB, then isSB = TRUE, otherwise FALSE
1777  // So, it is better, if one tests whether Gw is a GB
1778  // in ideals.cc:
1779  // idLift (ideal mod, ideal submod,ideal * rest, BOOLEAN goodShape,
1780  //           BOOLEAN isSB,BOOLEAN divide,matrix * unit)
1781
1782  // Let be Mtmp = {m1,...,ms}, where mi=sum hij.in_gj, for all i=1,...,s
1783  // We compute F = {f1,...,fs}, where fi=sum hij.gj
1784  int i, j, nM = IDELEMS(Mtmp);
1785  ideal idpol, idLG;
1786  ideal F = idInit(nM, 1);
1787
1788  for(i=0; i<nM; i++)
1789  {
1790     idpol = idVec2Ideal(Mtmp->m[i]);
1791     idLG = MidMult(idpol, G);
1792     idpol = NULL;
1793     F->m[i] = NULL;
1794     for(j=IDELEMS(idLG)-1; j>=0; j--)
1795     {
1796       F->m[i] = pAdd(F->m[i], idLG->m[j]);
1797       idLG->m[j]=NULL;
1798     }
1799     idDelete(&idLG);
1800  }
1801  idDelete(&Mtmp);
1802  return F;
1803}
1804
1805//unused
1806#if 0
1807static void checkidealCC(ideal G, char* Ch)
1808{
1809  int i,nmon=0,ntmp;
1810  int nG = IDELEMS(G);
1811  int n = nG-1;
1812  Print("\n//** Ideal %s besteht aus %d Polynomen mit ", Ch, nG);
1813
1814  for(i=0; i<nG; i++)
1815  {
1816    ntmp =  pLength(G->m[i]);
1817    nmon += ntmp;
1818
1819    if(i != n)
1820    {
1821      Print("%d, ", ntmp);
1822    }
1823    else
1824    {
1825      Print(" bzw. %d ", ntmp);
1826    }
1827  }
1828  PrintS(" Monomen.\n");
1829  Print("//** %s besitzt %d Monome.", Ch, nmon);
1830  PrintLn();
1831}
1832#endif
1833
1834//unused
1835#if 0
1836static void HeadidString(ideal L, char* st)
1837{
1838  int i, nL = IDELEMS(L)-1;
1839
1840  Print("//  The head terms of the ideal %s = ", st);
1841  for(i=0; i<nL; i++)
1842  {
1843    Print(" %s, ", pString(pHead(L->m[i])));
1844  }
1845  Print(" %s;\n", pString(pHead(L->m[nL])));
1846}
1847#endif
1848
1849static inline int MivComp(intvec* iva, intvec* ivb)
1850{
1851  assume(iva->length() == ivb->length());
1852  int i;
1853  for(i=iva->length()-1; i>=0; i--)
1854  {
1855    if((*iva)[i] - (*ivb)[i] != 0)
1856    {
1857      return 0;
1858    }
1859  }
1860  return 1;
1861}
1862
1863/**********************************************
1864 * Look for the smallest absolut value in vec *
1865 **********************************************/
1866static int MivAbsMax(intvec* vec)
1867{
1868  int i,k;
1869  if((*vec)[0] < 0)
1870  {
1871    k = -(*vec)[0];
1872  }
1873  else
1874  {
1875    k = (*vec)[0];
1876  }
1877  for(i=1; i < (vec->length()); i++)
1878  {
1879    if((*vec)[i] < 0)
1880    {
1881      if(-(*vec)[i] > k)
1882      {
1883        k = -(*vec)[i];
1884      }
1885    }
1886    else
1887    {
1888      if((*vec)[i] > k)
1889      {
1890        k = (*vec)[i];
1891      }
1892    }
1893  }
1894  return k;
1895}
1896
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
1917/**********************************************************************
1918 * Compute a next weight vector between curr_weight and target_weight *
1919 * with respect to an ideal <G>.                                      *
1920**********************************************************************/
1921static intvec* MwalkNextWeightCC(intvec* curr_weight, intvec* target_weight,
1922                                 ideal G)
1923{
1924  BOOLEAN nError = Overflow_Error;
1925  Overflow_Error = FALSE;
1926
1927  assume(currRing != NULL && curr_weight != NULL &&
1928         target_weight != NULL && G != NULL);
1929
1930  int nRing = currRing->N;
1931  int checkRed, j, nG = IDELEMS(G);
1932  intvec* ivtemp;
1933
1934  mpz_t t_zaehler, t_nenner;
1935  mpz_init(t_zaehler);
1936  mpz_init(t_nenner);
1937
1938  mpz_t s_zaehler, s_nenner, temp, MwWd;
1939  mpz_init(s_zaehler);
1940  mpz_init(s_nenner);
1941  mpz_init(temp);
1942  mpz_init(MwWd);
1943
1944  mpz_t sing_int;
1945  mpz_init(sing_int);
1946  mpz_set_si(sing_int,  2147483647);
1947
1948  mpz_t sing_int_half;
1949  mpz_init(sing_int_half);
1950  mpz_set_si(sing_int_half,  3*(1073741824/2));
1951
1952  mpz_t deg_w0_p1, deg_d0_p1;
1953  mpz_init(deg_w0_p1);
1954  mpz_init(deg_d0_p1);
1955
1956  mpz_t sztn, sntz;
1957  mpz_init(sztn);
1958  mpz_init(sntz);
1959
1960  mpz_t t_null;
1961  mpz_init(t_null);
1962
1963  mpz_t ggt;
1964  mpz_init(ggt);
1965
1966  mpz_t dcw;
1967  mpz_init(dcw);
1968
1969  int gcd_tmp;
1970  intvec* diff_weight = MivSub(target_weight, curr_weight);
1971
1972  intvec* diff_weight1 = MivSub(target_weight, curr_weight);
1973  poly g;
1974
1975  for (j=0; j<nG; j++)
1976  {
1977    g = G->m[j];
1978    if (g != NULL)
1979    {
1980      ivtemp = MExpPol(g);
1981      mpz_set_si(deg_w0_p1, MivDotProduct(ivtemp, curr_weight));
1982      mpz_set_si(deg_d0_p1, MivDotProduct(ivtemp, diff_weight));
1983      delete ivtemp;
1984
1985      pIter(g);
1986      while (g != NULL)
1987      {
1988        ivtemp = MExpPol(g);
1989        mpz_set_si(MwWd, MivDotProduct(ivtemp, curr_weight));
1990        mpz_sub(s_zaehler, deg_w0_p1, MwWd);
1991        if(mpz_cmp(s_zaehler, t_null) != 0)
1992        {
1993          mpz_set_si(MwWd, MivDotProduct(ivtemp, diff_weight));
1994          mpz_sub(s_nenner, MwWd, deg_d0_p1);
1995          // check for 0 < s <= 1
1996          if( (mpz_cmp(s_zaehler,t_null) > 0 &&
1997               mpz_cmp(s_nenner, s_zaehler)>=0) ||
1998              (mpz_cmp(s_zaehler, t_null) < 0 &&
1999               mpz_cmp(s_nenner, s_zaehler)<=0))
2000          {
2001            // make both positive
2002            if (mpz_cmp(s_zaehler, t_null) < 0)
2003            {
2004              mpz_neg(s_zaehler, s_zaehler);
2005              mpz_neg(s_nenner, s_nenner);
2006            }
2007
2008            //compute a simple fraction of s
2009            cancel(s_zaehler, s_nenner);
2010
2011            if(mpz_cmp(t_nenner, t_null) != 0)
2012            {
2013              mpz_mul(sztn, s_zaehler, t_nenner);
2014              mpz_mul(sntz, s_nenner, t_zaehler);
2015
2016              if(mpz_cmp(sztn,sntz) < 0)
2017              {
2018                mpz_add(t_nenner, t_null, s_nenner);
2019                mpz_add(t_zaehler,t_null, s_zaehler);
2020              }
2021            }
2022            else
2023            {
2024              mpz_add(t_nenner, t_null, s_nenner);
2025              mpz_add(t_zaehler,t_null, s_zaehler);
2026            }
2027          }
2028        }
2029        pIter(g);
2030        delete ivtemp;
2031      }
2032    }
2033  }
2034  //Print("\n// Alloc Size = %d \n", nRing*sizeof(mpz_t));
2035  mpz_t *vec=(mpz_t*)omAlloc(nRing*sizeof(mpz_t));
2036
2037
2038  // there is no 0<t<1 and define the next weight vector that is equal
2039  // to the current weight vector
2040  if(mpz_cmp(t_nenner, t_null) == 0)
2041  {
2042    #ifndef SING_NDEBUG
2043    Print("\n//MwalkNextWeightCC: t_nenner ist Null!");
2044    #endif
2045    delete diff_weight;
2046    diff_weight = ivCopy(curr_weight);//take memory
2047    goto FINISH;
2048  }
2049
2050  // define the target vector as the next weight vector, if t = 1
2051  if(mpz_cmp_si(t_nenner, 1)==0 && mpz_cmp_si(t_zaehler,1)==0)
2052  {
2053    delete diff_weight;
2054    diff_weight = ivCopy(target_weight); //this takes memory
2055    goto FINISH;
2056  }
2057
2058   checkRed = 0;
2059
2060  SIMPLIFY_GCD:
2061
2062  // simplify the vectors curr_weight and diff_weight (C-int)
2063  gcd_tmp = (*curr_weight)[0];
2064
2065  for (j=1; j<nRing; j++)
2066  {
2067    gcd_tmp = gcd(gcd_tmp, (*curr_weight)[j]);
2068    if(gcd_tmp == 1)
2069    {
2070      break;
2071    }
2072  }
2073  if(gcd_tmp != 1)
2074  {
2075    for (j=0; j<nRing; j++)
2076    {
2077      gcd_tmp = gcd(gcd_tmp, (*diff_weight)[j]);
2078      if(gcd_tmp == 1)
2079      {
2080        break;
2081      }
2082    }
2083  }
2084  if(gcd_tmp != 1)
2085  {
2086    for (j=0; j<nRing; j++)
2087    {
2088      (*curr_weight)[j] =  (*curr_weight)[j]/gcd_tmp;
2089      (*diff_weight)[j] =  (*diff_weight)[j]/gcd_tmp;
2090    }
2091  }
2092  if(checkRed > 0)
2093  {
2094    for (j=0; j<nRing; j++)
2095    {
2096      mpz_set_si(vec[j], (*diff_weight)[j]);
2097    }
2098    goto TEST_OVERFLOW;
2099  }
2100
2101#ifdef  NEXT_VECTORS_CC
2102  Print("\n// gcd of the weight vectors (current and target) = %d", gcd_tmp);
2103  ivString(curr_weight, "new cw");
2104  ivString(diff_weight, "new dw");
2105
2106  PrintS("\n// t_zaehler: ");  mpz_out_str( stdout, 10, t_zaehler);
2107  PrintS(", t_nenner: ");  mpz_out_str( stdout, 10, t_nenner);
2108#endif
2109
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
2116  for (j=0; j<nRing; j++)
2117  {
2118    mpz_set_si(dcw, (*curr_weight)[j]);
2119    mpz_mul(s_nenner, t_nenner, dcw);
2120
2121    if( (*diff_weight)[j]>0)
2122    {
2123      mpz_mul_ui(s_zaehler, t_zaehler, (*diff_weight)[j]);
2124    }
2125    else
2126    {
2127      mpz_mul_ui(s_zaehler, t_zaehler, -(*diff_weight)[j]);
2128      mpz_neg(s_zaehler, s_zaehler);
2129    }
2130    mpz_add(sntz, s_nenner, s_zaehler);
2131    mpz_init_set(vec[j], sntz);
2132
2133#ifdef NEXT_VECTORS_CC
2134    Print("\n//   j = %d ==> ", j);
2135    PrintS("(");
2136    mpz_out_str( stdout, 10, t_nenner);
2137    Print(" * %d)", (*curr_weight)[j]);
2138    Print(" + ("); mpz_out_str( stdout, 10, t_zaehler);
2139    Print(" * %d) =  ",  (*diff_weight)[j]);
2140    mpz_out_str( stdout, 10, s_nenner);
2141    PrintS(" + ");
2142    mpz_out_str( stdout, 10, s_zaehler);
2143    PrintS(" = "); mpz_out_str( stdout, 10, sntz);
2144    Print(" ==> vector[%d]: ", j); mpz_out_str(stdout, 10, vec[j]);
2145#endif
2146
2147    if(j==0)
2148    {
2149      mpz_set(ggt, sntz);
2150    }
2151    else
2152    {
2153      if(mpz_cmp_si(ggt,1) != 0)
2154      {
2155        mpz_gcd(ggt, ggt, sntz);
2156      }
2157    }
2158  }
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  }
2167#ifdef  NEXT_VECTORS_CC
2168  PrintS("\n// gcd of elements of the vector: ");
2169  mpz_out_str( stdout, 10, ggt);
2170#endif
2171
2172  for(j=0; j<nRing; j++)
2173    {
2174    if(mpz_cmp(vec[j], sing_int_half) >= 0)
2175      {
2176      goto REDUCTION;
2177      }
2178    }
2179  checkRed = 1;
2180  for (j=0; j<nRing; j++)
2181    {
2182      (*diff_weight)[j] = mpz_get_si(vec[j]);
2183    }
2184  goto SIMPLIFY_GCD;
2185
2186  REDUCTION:
2187  checkRed = 1;
2188  for (j=0; j<nRing; j++)
2189  {
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;
2219  }
2220
2221 TEST_OVERFLOW:
2222
2223  for (j=0; j<nRing; j++)
2224  {
2225    if(mpz_cmp(vec[j], sing_int)>=0)
2226    {
2227      if(Overflow_Error == FALSE)
2228      {
2229        Overflow_Error = TRUE;
2230        PrintS("\n// ** OVERFLOW in \"MwalkNextWeightCC\": ");
2231        mpz_out_str( stdout, 10, vec[j]);
2232        PrintS(" is greater than 2147483647 (max. integer representation)\n");
2233        //Print("//  So vector[%d] := %d is wrong!!\n",j+1, vec[j]);// vec[j] is mpz_t
2234      }
2235    }
2236  }
2237
2238 FINISH:
2239   delete diff_weight1;
2240   mpz_clear(t_zaehler);
2241   mpz_clear(t_nenner);
2242   mpz_clear(s_zaehler);
2243   mpz_clear(s_nenner);
2244   mpz_clear(sntz);
2245   mpz_clear(sztn);
2246   mpz_clear(temp);
2247   mpz_clear(MwWd);
2248   mpz_clear(deg_w0_p1);
2249   mpz_clear(deg_d0_p1);
2250   mpz_clear(ggt);
2251   omFree(vec);
2252   mpz_clear(sing_int_half);
2253   mpz_clear(sing_int);
2254   mpz_clear(dcw);
2255   mpz_clear(t_null);
2256
2257  if(Overflow_Error == FALSE)
2258  {
2259    Overflow_Error = nError;
2260  }
2261  rComplete(currRing);
2262  for(j=0; j<IDELEMS(G); j++)
2263  {
2264    poly p=G->m[j];
2265    while(p!=NULL)
2266    {
2267      p_Setm(p,currRing);
2268      pIter(p);
2269    }
2270  }
2271return diff_weight;
2272}
2273
2274/**********************************************************************
2275* Compute an intermediate weight vector from iva to ivb w.r.t.        *
2276* the reduced Groebner basis G.                                       *
2277* Return NULL, if it is equal to iva or iva = avb.                    *
2278**********************************************************************/
2279intvec* MkInterRedNextWeight(intvec* iva, intvec* ivb, ideal G)
2280{
2281  intvec* tmp = new intvec(iva->length());
2282  intvec* result;
2283
2284  if(G == NULL)
2285  {
2286    return tmp;
2287  }
2288  if(MivComp(iva, ivb) == 1)
2289  {
2290    return tmp;
2291  }
2292  result = MwalkNextWeightCC(iva, ivb, G);
2293
2294  if(MivComp(result, iva) == 1)
2295  {
2296    delete result;
2297    return tmp;
2298  }
2299
2300  delete tmp;
2301  return result;
2302}
2303
2304/********************************************************************
2305 * define and execute a new ring which order is (a(vb),a(va),lp,C)  *
2306 * ******************************************************************/
2307static void VMrHomogeneous(intvec* va, intvec* vb)
2308{
2309
2310  if ((currRing->ppNoether)!=NULL)
2311  {
2312    pDelete(&(currRing->ppNoether));
2313  }
2314  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
2315      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
2316  {
2317    sLastPrinted.CleanUp();
2318  }
2319
2320  ring r = (ring) omAlloc0Bin(sip_sring_bin);
2321  int i, nv = currRing->N;
2322
2323  r->cf  = currRing->cf;
2324  r->N   = currRing->N;
2325  int nb = 4;
2326
2327
2328  //names
2329  char* Q; // In order to avoid the corrupted memory, do not change.
2330  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2331  for(i=0; i<nv; i++)
2332  {
2333    Q = currRing->names[i];
2334    r->names[i]  = omStrDup(Q);
2335  }
2336
2337  //weights: entries for 3 blocks: NULL Made:???
2338  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2339  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
2340  r->wvhdl[1] = (int*) omAlloc((nv-1)*sizeof(int));
2341
2342  for(i=0; i<nv-1; i++)
2343  {
2344    r->wvhdl[1][i] = (*vb)[i];
2345    r->wvhdl[0][i] = (*va)[i];
2346  }
2347  r->wvhdl[0][nv] = (*va)[nv];
2348
2349  // order: (1..1),a,lp,C
2350  r->order = (int *) omAlloc(nb * sizeof(int *));
2351  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2352  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2353
2354  // ringorder a for the first block: var 1..nv
2355  r->order[0]  = ringorder_a;
2356  r->block0[0] = 1;
2357  r->block1[0] = nv;
2358
2359 // ringorder a for the second block: var 2..nv
2360  r->order[1]  = ringorder_a;
2361  r->block0[1] = 2;
2362  r->block1[1] = nv;
2363
2364  // ringorder lp for the third block: var 2..nv
2365  r->order[2]  = ringorder_lp;
2366  r->block0[2] = 2;
2367  r->block1[2] = nv;
2368
2369  // ringorder C for the 4th block
2370  // it is very important within "idLift",
2371  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
2372  // therefore, nb  must be (nBlocks(currRing)  + 1)
2373  r->order[3]  = ringorder_C;
2374
2375  // polynomial ring
2376  r->OrdSgn    = 1;
2377
2378  // complete ring intializations
2379
2380  rComplete(r);
2381
2382  rChangeCurrRing(r);
2383}
2384
2385
2386/**************************************************************
2387 * define and execute a new ring which order is (a(va),lp,C)  *
2388 * ************************************************************/
2389static ring VMrDefault(intvec* va)
2390{
2391
2392  if ((currRing->ppNoether)!=NULL)
2393  {
2394    pDelete(&(currRing->ppNoether));
2395  }
2396  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
2397      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
2398  {
2399    sLastPrinted.CleanUp();
2400  }
2401
2402  ring r = (ring) omAlloc0Bin(sip_sring_bin);
2403  int i, nv = currRing->N;
2404
2405  r->cf  = currRing->cf;
2406  r->N   = currRing->N;
2407
2408  int nb = 4;
2409
2410  //names
2411  char* Q; // In order to avoid the corrupted memory, do not change.
2412  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2413  for(i=0; i<nv; i++)
2414  {
2415    Q = currRing->names[i];
2416    r->names[i]  = omStrDup(Q);
2417  }
2418
2419  /*weights: entries for 3 blocks: NULL Made:???*/
2420  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2421  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
2422  for(i=0; i<nv; i++)
2423    r->wvhdl[0][i] = (*va)[i];
2424
2425  /* order: a,lp,C,0 */
2426  r->order = (int *) omAlloc(nb * sizeof(int *));
2427  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2428  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2429
2430  // ringorder a for the first block: var 1..nv
2431  r->order[0]  = ringorder_a;
2432  r->block0[0] = 1;
2433  r->block1[0] = nv;
2434
2435  // ringorder lp for the second block: var 1..nv
2436  r->order[1]  = ringorder_lp;
2437  r->block0[1] = 1;
2438  r->block1[1] = nv;
2439
2440  // ringorder C for the third block
2441  // it is very important within "idLift",
2442  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
2443  // therefore, nb  must be (nBlocks(currRing)  + 1)
2444  r->order[2]  = ringorder_C;
2445
2446  // the last block: everything is 0
2447  r->order[3]  = 0;
2448
2449  // polynomial ring
2450  r->OrdSgn    = 1;
2451
2452  // complete ring intializations
2453
2454  rComplete(r);
2455  return r;
2456  //rChangeCurrRing(r);
2457}
2458
2459/****************************************************************
2460 * define and execute a new ring with ordering (a(va),Wp(vb),C) *
2461 * **************************************************************/
2462static ring VMrRefine(intvec* va, intvec* vb)
2463{
2464
2465  if ((currRing->ppNoether)!=NULL)
2466  {
2467    pDelete(&(currRing->ppNoether));
2468  }
2469  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
2470      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
2471  {
2472    sLastPrinted.CleanUp();
2473  }
2474
2475  ring r = (ring) omAlloc0Bin(sip_sring_bin);
2476  int i, nv = currRing->N;
2477
2478  r->cf  = currRing->cf;
2479  r->N   = currRing->N;
2480  //int nb = nBlocks(currRing)  + 1;
2481  int nb = 4;
2482
2483  //names
2484  char* Q; // In order to avoid the corrupted memory, do not change.
2485  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2486  for(i=0; i<nv; i++)
2487  {
2488    Q = currRing->names[i];
2489    r->names[i]  = omStrDup(Q);
2490  }
2491
2492  //weights: entries for 3 blocks: NULL Made:???
2493  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2494  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
2495  r->wvhdl[1] = (int*) omAlloc(nv*sizeof(int));
2496
2497  for(i=0; i<nv; i++)
2498  {
2499    r->wvhdl[0][i] = (*va)[i];
2500    r->wvhdl[1][i] = (*vb)[i];
2501  }
2502  r->wvhdl[2]=NULL;
2503  r->wvhdl[3]=NULL;
2504
2505  // order: (1..1),a,lp,C
2506  r->order = (int *) omAlloc(nb * sizeof(int *));
2507  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2508  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2509
2510  // ringorder a for the first block: var 1..nv
2511  r->order[0]  = ringorder_a;
2512  r->block0[0] = 1;
2513  r->block1[0] = nv;
2514
2515 // ringorder Wp for the second block: var 1..nv
2516  r->order[1]  = ringorder_Wp;
2517  r->block0[1] = 1;
2518  r->block1[1] = nv;
2519
2520  // ringorder lp for the third block: var 1..nv
2521  r->order[2]  = ringorder_C;
2522  r->block0[2] = 1;
2523  r->block1[2] = nv;
2524
2525  // ringorder C for the 4th block
2526  // it is very important within "idLift",
2527  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
2528  // therefore, nb  must be (nBlocks(currRing)  + 1)
2529  r->order[3]  = 0;
2530
2531  // polynomial ring
2532  r->OrdSgn    = 1;
2533
2534  // complete ring intializations
2535 
2536  rComplete(r);
2537
2538  //rChangeCurrRing(r);
2539  return r;
2540}
2541
2542/*****************************************************
2543 * define and execute a new ring with ordering (M,C) *
2544 *****************************************************/
2545static ring VMatrDefault(intvec* va)
2546{
2547
2548  if ((currRing->ppNoether)!=NULL)
2549  {
2550    pDelete(&(currRing->ppNoether));
2551  }
2552  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
2553      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
2554  {
2555    sLastPrinted.CleanUp();
2556  }
2557
2558  ring r = (ring) omAlloc0Bin(sip_sring_bin);
2559  int i, nv = currRing->N;
2560
2561  r->cf  = currRing->cf;
2562  r->N   = currRing->N;
2563
2564  int nb = 4;
2565
2566  //names
2567  char* Q; // In order to avoid the corrupted memory, do not change.
2568  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2569  for(i=0; i<nv; i++)
2570  {
2571    Q = currRing->names[i];
2572    r->names[i]  = omStrDup(Q);
2573  }
2574
2575  /*weights: entries for 3 blocks: NULL Made:???*/
2576  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2577  r->wvhdl[0] = (int*) omAlloc(nv*nv*sizeof(int));
2578  r->wvhdl[1] =NULL; // (int*) omAlloc(nv*sizeof(int));
2579  r->wvhdl[2]=NULL;
2580  r->wvhdl[3]=NULL;
2581  for(i=0; i<nv*nv; i++)
2582    r->wvhdl[0][i] = (*va)[i];
2583
2584  /* order: a,lp,C,0 */
2585  r->order = (int *) omAlloc(nb * sizeof(int *));
2586  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2587  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2588
2589  // ringorder a for the first block: var 1..nv
2590  r->order[0]  = ringorder_M;
2591  r->block0[0] = 1;
2592  r->block1[0] = nv;
2593
2594  // ringorder C for the second block
2595  r->order[1]  = ringorder_C;
2596  r->block0[1] = 1;
2597  r->block1[1] = nv;
2598
2599
2600// ringorder C for the third block: var 1..nv
2601  r->order[2]  = ringorder_C;
2602  r->block0[2] = 1;
2603  r->block1[2] = nv;
2604
2605
2606  // the last block: everything is 0
2607  r->order[3]  = 0;
2608
2609  // polynomial ring
2610  r->OrdSgn    = 1;
2611
2612  // complete ring intializations
2613
2614  rComplete(r);
2615
2616  //rChangeCurrRing(r);
2617  return r;
2618}
2619
2620/***********************************************************
2621 * define and execute a new ring with ordering (a(vb),M,C) *
2622 ***********************************************************/
2623static ring VMatrRefine(intvec* va, intvec* vb)
2624{
2625
2626  if ((currRing->ppNoether)!=NULL)
2627  {
2628    pDelete(&(currRing->ppNoether));
2629  }
2630  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
2631      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
2632  {
2633    sLastPrinted.CleanUp();
2634  }
2635
2636  ring r = (ring) omAlloc0Bin(sip_sring_bin);
2637  int i, nv = currRing->N;
2638  int nvs = nv*nv;
2639  r->cf  = currRing->cf;
2640  r->N   = currRing->N;
2641
2642  int nb = 4;
2643
2644  //names
2645  char* Q; // In order to avoid the corrupted memory, do not change.
2646  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2647  for(i=0; i<nv; i++)
2648  {
2649    Q = currRing->names[i];
2650    r->names[i]  = omStrDup(Q);
2651  }
2652
2653  /*weights: entries for 3 blocks: NULL Made:???*/
2654  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2655  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
2656  r->wvhdl[1] = (int*) omAlloc(nvs*sizeof(int));
2657  r->wvhdl[2]=NULL;
2658  r->wvhdl[3]=NULL;
2659  for(i=0; i<nvs; i++)
2660  {
2661    r->wvhdl[1][i] = (*va)[i];
2662  }
2663  for(i=0; i<nv; i++)
2664  {
2665    r->wvhdl[0][i] = (*vb)[i];
2666  }
2667  /* order: a,lp,C,0 */
2668  r->order = (int *) omAlloc(nb * sizeof(int *));
2669  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2670  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2671
2672  // ringorder a for the first block: var 1..nv
2673  r->order[0]  = ringorder_a;
2674  r->block0[0] = 1;
2675  r->block1[0] = nv;
2676
2677  // ringorder M for the second block: var 1..nv
2678  r->order[1]  = ringorder_M;
2679  r->block0[1] = 1;
2680  r->block1[1] = nv;
2681
2682  // ringorder C for the third block: var 1..nv
2683  r->order[2]  = ringorder_C;
2684  r->block0[2] = 1;
2685  r->block1[2] = nv;
2686
2687
2688  // the last block: everything is 0
2689  r->order[3]  = 0;
2690
2691  // polynomial ring
2692  r->OrdSgn    = 1;
2693
2694  // complete ring intializations
2695
2696  rComplete(r);
2697
2698  //rChangeCurrRing(r);
2699  return r;
2700}
2701
2702/**********************************************************************
2703* define and execute a new ring which order is  a lexicographic order *
2704***********************************************************************/
2705static void VMrDefaultlp(void)
2706{
2707
2708  if ((currRing->ppNoether)!=NULL)
2709  {
2710    pDelete(&(currRing->ppNoether));
2711  }
2712  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
2713      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
2714
2715  {
2716    sLastPrinted.CleanUp();
2717  }
2718
2719  ring r = (ring) omAlloc0Bin(sip_sring_bin);
2720  int i, nv = currRing->N;
2721
2722  r->cf  = currRing->cf;
2723  r->N   = currRing->N;
2724  int nb = rBlocks(currRing) + 1;
2725
2726  // names
2727  char* Q; // to avoid the corrupted memory, do not change!!
2728  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2729  for(i=0; i<nv; i++)
2730  {
2731    Q = currRing->names[i];
2732    r->names[i]  = omStrDup(Q);
2733  }
2734
2735  /*weights: entries for 3 blocks: NULL Made:???*/
2736
2737  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2738
2739  /* order: lp,C,0 */
2740  r->order = (int *) omAlloc(nb * sizeof(int *));
2741  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2742  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2743
2744  /* ringorder lp for the first block: var 1..nv */
2745  r->order[0]  = ringorder_lp;
2746  r->block0[0] = 1;
2747  r->block1[0] = nv;
2748
2749  /* ringorder C for the second block */
2750  r->order[1]  = ringorder_C;
2751
2752  /* the last block: everything is 0 */
2753  r->order[2]  = 0;
2754
2755  /*polynomial ring*/
2756  r->OrdSgn    = 1;
2757
2758  /* complete ring intializations */
2759
2760  rComplete(r);
2761
2762  rChangeCurrRing(r);
2763}
2764
2765/***************************************************
2766* define a ring with parameters und change to it   *
2767* DefRingPar and DefRingParlp corrupt still memory *
2768****************************************************/
2769static void DefRingPar(intvec* va)
2770{
2771  int i, nv = currRing->N;
2772  int nb = rBlocks(currRing) + 1;
2773
2774  ring res=(ring)omAllocBin(sip_sring_bin);
2775
2776  memcpy(res,currRing,sizeof(ip_sring));
2777
2778  res->VarOffset = NULL;
2779  res->ref=0;
2780
2781  res->cf = currRing->cf; currRing->cf->ref++;
2782
2783
2784  /*weights: entries for 3 blocks: NULL Made:???*/
2785  res->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2786  res->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
2787  for(i=0; i<nv; i++)
2788    res->wvhdl[0][i] = (*va)[i];
2789
2790  /* order: a,lp,C,0 */
2791
2792  res->order = (int *) omAlloc(nb * sizeof(int *));
2793  res->block0 = (int *)omAlloc0(nb * sizeof(int *));
2794  res->block1 = (int *)omAlloc0(nb * sizeof(int *));
2795
2796  // ringorder a for the first block: var 1..nv
2797  res->order[0]  = ringorder_a;
2798  res->block0[0] = 1;
2799  res->block1[0] = nv;
2800
2801  // ringorder lp for the second block: var 1..nv
2802  res->order[1]  = ringorder_lp;
2803  res->block0[1] = 1;
2804  res->block1[1] = nv;
2805
2806  // ringorder C for the third block
2807  // it is very important within "idLift",
2808  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
2809  // therefore, nb  must be (nBlocks(currRing)  + 1)
2810  res->order[2]  = ringorder_C;
2811
2812  // the last block: everything is 0
2813  res->order[3]  = 0;
2814
2815  // polynomial ring
2816  res->OrdSgn    = 1;
2817
2818
2819  res->names   = (char **)omAlloc0(nv * sizeof(char_ptr));
2820  for (i=nv-1; i>=0; i--)
2821  {
2822    res->names[i] = omStrDup(currRing->names[i]);
2823  }
2824  // complete ring intializations
2825  rComplete(res);
2826
2827  // clean up history
2828  if (sLastPrinted.RingDependend())
2829  {
2830    sLastPrinted.CleanUp();
2831  }
2832
2833
2834  // execute the created ring
2835  rChangeCurrRing(res);
2836}
2837
2838
2839static void DefRingParlp(void)
2840{
2841  int i, nv = currRing->N;
2842
2843  ring r=(ring)omAllocBin(sip_sring_bin);
2844
2845  memcpy(r,currRing,sizeof(ip_sring));
2846
2847  r->VarOffset = NULL;
2848  r->ref=0;
2849
2850  r->cf = currRing->cf; currRing->cf->ref++;
2851
2852
2853  r->cf  = currRing->cf;
2854  r->N   = currRing->N;
2855  int nb = rBlocks(currRing) + 1;
2856
2857  // names
2858  char* Q;
2859  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2860  for(i=nv-1; i>=0; i--)
2861  {
2862    Q = currRing->names[i];
2863    r->names[i]  = omStrDup(Q);
2864  }
2865
2866  /*weights: entries for 3 blocks: NULL Made:???*/
2867
2868  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2869
2870  /* order: lp,C,0 */
2871  r->order = (int *) omAlloc(nb * sizeof(int *));
2872  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2873  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2874
2875  /* ringorder lp for the first block: var 1..nv */
2876  r->order[0]  = ringorder_lp;
2877  r->block0[0] = 1;
2878  r->block1[0] = nv;
2879
2880  /* ringorder C for the second block */
2881  r->order[1]  = ringorder_C;
2882
2883  /* the last block: everything is 0 */
2884  r->order[2]  = 0;
2885
2886  /*polynomial ring*/
2887  r->OrdSgn    = 1;
2888
2889
2890//   if (rParameter(currRing)!=NULL)
2891//   {
2892//     r->cf->extRing->qideal->m[0]=p_Copy(currRing->cf->extRing->qideal->m[0], currRing->cf->extRing);
2893//     int l=rPar(currRing);
2894//     r->cf->extRing->names=(char **)omAlloc(l*sizeof(char_ptr));
2895//
2896//     for(i=l-1;i>=0;i--)
2897//     {
2898//       rParameter(r)[i]=omStrDup(rParameter(currRing)[i]);
2899//     }
2900//   }
2901
2902  // complete ring intializations
2903
2904  rComplete(r);
2905
2906  // clean up history
2907  if (sLastPrinted.RingDependend())
2908  {
2909    sLastPrinted.CleanUp();
2910  }
2911
2912  // execute the created ring
2913  rChangeCurrRing(r);
2914}
2915
2916/*************************************************************
2917 * check whether one or more components of a vector are zero *
2918 *************************************************************/
2919static int isNolVector(intvec* hilb)
2920{
2921  int i;
2922  for(i=hilb->length()-1; i>=0; i--)
2923  {
2924    if((* hilb)[i]==0)
2925    {
2926      return 1;
2927    }
2928  }
2929  return 0;
2930}
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}
3023
3024/******************************  Februar 2002  ****************************
3025 * G is a Groebner basis w.r.t. (a(curr_weight),lp) and                   *
3026 * we compute a GB of <G> w.r.t. the lex. order by the perturbation walk  *
3027 * its perturbation degree is tp_deg                                      *
3028 * We call the following subfunction LastGB, if                           *
3029 * the computed intermediate weight vector or                             *
3030 * if the perturbed target weight vector does NOT lie n the correct cone  *
3031 **************************************************************************/
3032
3033static ideal LastGB(ideal G, intvec* curr_weight,int tp_deg)
3034{
3035  BOOLEAN nError = Overflow_Error;
3036  Overflow_Error = FALSE;
3037
3038  int i, nV = currRing->N;
3039  int nwalk=0, endwalks=0, nnwinC=1;
3040  int nlast = 0;
3041  ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG;
3042  ring newRing, oldRing, TargetRing;
3043  intvec* iv_M_lp;
3044  intvec* target_weight;
3045  intvec* iv_lp = Mivlp(nV); //define (1,0,...,0)
3046  intvec* pert_target_vector;
3047  intvec* ivNull = new intvec(nV);
3048  intvec* extra_curr_weight = new intvec(nV);
3049  intvec* next_weight;
3050
3051#ifndef  BUCHBERGER_ALG
3052  intvec* hilb_func;
3053#endif
3054
3055  // to avoid (1,0,...,0) as the target vector
3056  intvec* last_omega = new intvec(nV);
3057  for(i=nV-1; i>0; i--)
3058  {
3059    (*last_omega)[i] = 1;
3060  }
3061  (*last_omega)[0] = 10000;
3062
3063  ring EXXRing = currRing;
3064
3065  // compute a pertubed weight vector of the target weight vector
3066  if(tp_deg > 1 && tp_deg <= nV)
3067  {
3068    //..25.03.03 VMrDefaultlp();//    VMrDefault(target_weight);
3069    if (rParameter (currRing) != NULL)
3070    {
3071      DefRingParlp();
3072    }
3073    else
3074    {
3075      VMrDefaultlp();
3076    }
3077    TargetRing = currRing;
3078    ssG = idrMoveR(G,EXXRing,currRing);
3079    iv_M_lp = MivMatrixOrderlp(nV);
3080    //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
3081    target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
3082    delete iv_M_lp;
3083    pert_target_vector = target_weight;
3084
3085    rChangeCurrRing(EXXRing);
3086    G = idrMoveR(ssG, TargetRing,currRing);
3087  }
3088  else
3089  {
3090    target_weight = Mivlp(nV);
3091  }
3092  //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
3093
3094  while(1)
3095  {
3096    nwalk++;
3097    nstep++;
3098    to=clock();
3099   // compute a next weight vector
3100    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
3101    xtnw=xtnw+clock()-to;
3102
3103#ifdef PRINT_VECTORS
3104    MivString(curr_weight, target_weight, next_weight);
3105#endif
3106
3107    if(Overflow_Error == TRUE)
3108    {
3109      newRing = currRing;
3110      nnwinC = 0;
3111      if(tp_deg == 1)
3112      {
3113        nlast = 1;
3114      }
3115      delete next_weight;
3116
3117      //idElements(G, "G");
3118      //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
3119
3120      break;
3121    }
3122
3123    if(MivComp(next_weight, ivNull) == 1)
3124    {
3125      //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
3126      newRing = currRing;
3127      delete next_weight;
3128      break;
3129    }
3130
3131    if(MivComp(next_weight, target_weight) == 1)
3132      endwalks = 1;
3133
3134    for(i=nV-1; i>=0; i--)
3135    {
3136      (*extra_curr_weight)[i] = (*curr_weight)[i];
3137    }
3138    /* 06.11.01 NOT Changed */
3139    for(i=nV-1; i>=0; i--)
3140    {
3141      (*curr_weight)[i] = (*next_weight)[i];
3142    }
3143    oldRing = currRing;
3144    to=clock();
3145    // compute an initial form ideal of <G> w.r.t. "curr_vector"
3146    Gomega = MwalkInitialForm(G, curr_weight);
3147    xtif=xtif+clock()-to;
3148
3149#ifdef ENDWALKS
3150    if(endwalks == 1)
3151    {
3152      Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
3153      idElements(Gomega, "Gw");
3154      headidString(Gomega, "Gw");
3155    }
3156#endif
3157
3158#ifndef  BUCHBERGER_ALG
3159    if(isNolVector(curr_weight) == 0)
3160    {
3161      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
3162    }
3163    else
3164    {
3165      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
3166    }
3167#endif // BUCHBERGER_ALG
3168
3169    /* define a new ring that its ordering is "(a(curr_weight),lp) */
3170    //..25.03.03 VMrDefault(curr_weight);
3171    if (rParameter (currRing) != NULL)
3172    {
3173      DefRingPar(curr_weight);
3174    }
3175    else
3176    {
3177      rChangeCurrRing(VMrDefault(curr_weight));
3178    }
3179    newRing = currRing;
3180    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
3181
3182    to=clock();
3183    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
3184#ifdef  BUCHBERGER_ALG
3185    M = MstdhomCC(Gomega1);
3186#else
3187    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
3188    delete hilb_func;
3189#endif // BUCHBERGER_ALG
3190    xtstd=xtstd+clock()-to;
3191    /* change the ring to oldRing */
3192    rChangeCurrRing(oldRing);
3193    M1 =  idrMoveR(M, newRing,currRing);
3194    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
3195
3196    to=clock();
3197    /* compute a reduced Groebner basis of <G> w.r.t. "newRing" */
3198    F = MLifttwoIdeal(Gomega2, M1, G);
3199    xtlift=xtlift+clock()-to;
3200
3201    idDelete(&M1);
3202    idDelete(&G);
3203
3204    /* change the ring to newRing */
3205    rChangeCurrRing(newRing);
3206    F1 = idrMoveR(F, oldRing,currRing);
3207
3208    to=clock();
3209    /* reduce the Groebner basis <G> w.r.t. new ring */
3210    G = kInterRedCC(F1, NULL);
3211    xtred=xtred+clock()-to;
3212    idDelete(&F1);
3213
3214    if(endwalks == 1)
3215    {
3216      //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
3217      break;
3218    }
3219
3220    delete next_weight;
3221  }//while
3222
3223  delete ivNull;
3224
3225  if(tp_deg != 1)
3226  {
3227    //..25.03.03 VMrDefaultlp();//define and execute the ring "lp"
3228    if (rParameter (currRing) != NULL)
3229    {
3230      DefRingParlp();
3231    }
3232    else
3233    {
3234      VMrDefaultlp();
3235    }
3236    F1 = idrMoveR(G, newRing,currRing);
3237
3238    if(nnwinC == 0 || test_w_in_ConeCC(F1, pert_target_vector) != 1)
3239    {
3240      oldRing = currRing;
3241      rChangeCurrRing(newRing);
3242      G = idrMoveR(F1, oldRing,currRing);
3243      Print("\n// takes %d steps and calls the recursion of level %d:",
3244             nwalk, tp_deg-1);
3245
3246      F1 = LastGB(G,curr_weight, tp_deg-1);
3247    }
3248
3249    TargetRing = currRing;
3250    rChangeCurrRing(EXXRing);
3251    result = idrMoveR(F1, TargetRing,currRing);
3252  }
3253  else
3254  {
3255    if(nlast == 1)
3256    {
3257      //OMEGA_OVERFLOW_LASTGB:
3258      /*
3259      if(MivSame(curr_weight, iv_lp) == 1)
3260        if (rParameter(currRing) != NULL)
3261          DefRingParlp();
3262        else
3263          VMrDefaultlp();
3264      else
3265        if (rParameter(currRing) != NULL)
3266          DefRingPar(curr_weight);
3267        else
3268          VMrDefault(curr_weight);
3269      */
3270
3271        //..25.03.03 VMrDefaultlp();//define and execute the ring "lp"
3272        if (rParameter (currRing) != NULL)
3273        {
3274          DefRingParlp();
3275        }
3276        else
3277        {
3278          VMrDefaultlp();
3279        }
3280
3281      F1 = idrMoveR(G, newRing,currRing);
3282      //Print("\n// Apply \"std\" in ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
3283
3284      G = MstdCC(F1);
3285      idDelete(&F1);
3286      newRing = currRing;
3287    }
3288
3289    rChangeCurrRing(EXXRing);
3290    result = idrMoveR(G, newRing,currRing);
3291  }
3292  delete target_weight;
3293  delete last_omega;
3294  delete iv_lp;
3295
3296  if(Overflow_Error == FALSE)
3297  {
3298    Overflow_Error = nError;
3299  }
3300  return(result);
3301}
3302
3303/**********************************************************
3304 * check whether a polynomial of G has least 3 monomials  *
3305 **********************************************************/
3306static int lengthpoly(ideal G)
3307{
3308  int i;
3309  for(i=IDELEMS(G)-1; i>=0; i--)
3310  {
3311#if 0
3312    if(pLength(G->m[i])>2)
3313    {
3314      return 1;
3315    }
3316#else
3317    if((G->m[i]!=NULL) /* len >=0 */
3318       && (G->m[i]->next!=NULL) /* len >=1 */
3319       && (G->m[i]->next->next!=NULL) /* len >=2 */
3320       && (G->m[i]->next->next->next!=NULL) /* len >=3 */
3321      //&& (G->m[i]->next->next->next->next!=NULL) /* len >=4 */
3322       )
3323    {
3324    return 1;
3325    }
3326#endif
3327  }
3328  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;
3346}
3347
3348/*********************************************************
3349 * check whether a polynomial of G has least 2 monomials *
3350**********************************************************/
3351static int islengthpoly2(ideal G)
3352{
3353  int i;
3354  for(i=IDELEMS(G)-1; i>=0; i--)
3355  {
3356    if((G->m[i]!=NULL) /* len >=0 */
3357       && (G->m[i]->next!=NULL) /* len >=1 */
3358       && (G->m[i]->next->next!=NULL)) /* len >=2 */
3359    {
3360      return 1;
3361    }
3362  }
3363  return 0;
3364}
3365
3366
3367
3368/* Implementation of the improved Groebner walk algorithm which is written
3369   by Quoc-Nam Tran (2000).
3370   One perturbs the original target weight vector, only if
3371   the next intermediate weight vector is equal to the current target weight
3372   vector. This must be repeated until the wanted reduced Groebner basis
3373   to reach.
3374   If the numbers of variables is big enough, the representation of the origin
3375   weight vector may be very big. Therefore, it is possible the intermediate
3376   weight vector doesn't stay in the correct Groebner cone.
3377   In this case we have just a reduced Groebner basis of the given ideal
3378   with respect to another monomial order. Then we have to compute
3379   a wanted reduced Groebner basis of it with respect to the given order.
3380   At the following subroutine we use the improved Buchberger algorithm or
3381   the changed perturbation walk algorithm with a decrased degree.
3382 */
3383
3384/***************************************
3385 * return the initial term of an ideal *
3386 ***************************************/
3387static ideal idHeadCC(ideal h)
3388{
3389  int i, nH =IDELEMS(h);
3390
3391  ideal m = idInit(nH,h->rank);
3392
3393  for (i=nH-1;i>=0; i--)
3394  {
3395    if (h->m[i]!=NULL)
3396    {
3397      m->m[i]=pHead(h->m[i]);
3398    }
3399  }
3400  return m;
3401}
3402
3403/**********************************************
3404 * check whether two head-ideals are the same *
3405 **********************************************/
3406static inline int test_G_GB_walk(ideal H0, ideal H1)
3407{
3408  int i, nG = IDELEMS(H0);
3409
3410  if(nG != IDELEMS(H1))
3411  {
3412    return 0;
3413  }
3414  for(i=nG-1; i>=0; i--)
3415  {
3416#if 0
3417    poly t;
3418    if((t=pSub(pCopy(H0->m[i]), pCopy(H1->m[i]))) != NULL)
3419    {
3420      pDelete(&t);
3421      return 0;
3422    }
3423    pDelete(&t);
3424#else
3425    if(!pEqualPolys(H0->m[i],H1->m[i]))
3426    {
3427      return 0;
3428    }
3429#endif
3430  }
3431  return 1;
3432}
3433
3434//unused
3435/*****************************************************
3436 * find the maximal total degree of polynomials in G *
3437 *****************************************************/
3438#if 0
3439static int Trandegreebound(ideal G)
3440{
3441  int i, nG = IDELEMS(G);
3442  // int np=1;
3443  int nV = currRing->N;
3444  int degtmp, result = 0;
3445  intvec* ivUnit = Mivdp(nV);
3446
3447  for(i=nG-1; i>=0; i--)
3448  {
3449    // find the maximal total degree of the polynomial G[i]
3450    degtmp = MwalkWeightDegree(G->m[i], ivUnit);
3451    if(degtmp > result)
3452    {
3453      result = degtmp;
3454    }
3455  }
3456  delete ivUnit;
3457  return result;
3458}
3459#endif
3460
3461//unused
3462/************************************************************************
3463 * perturb the weight vector iva w.r.t. the ideal G.                    *
3464 *  the monomial order of the current ring is the w_1 weight lex. order *
3465 *  define w := d^(n-1)w_1+ d^(n-2)w_2, ...+ dw_(n-1)+ w_n              *
3466 *  where d := 1 + max{totdeg(g):g in G}*m, or                          *
3467 *  d := (2*maxdeg*maxdeg + (nV+1)*maxdeg)*m;                           *
3468 ************************************************************************/
3469#if 0
3470static intvec* TranPertVector(ideal G, intvec* iva)
3471{
3472  BOOLEAN nError = Overflow_Error;
3473  Overflow_Error = FALSE;
3474
3475  int i, j;
3476  // int nG = IDELEMS(G);
3477  int nV = currRing->N;
3478
3479  // define the sequence which expresses the current monomial ordering
3480  // w_1 = iva; w_2 = (1,0,..,0); w_n = (0,...,0,1,0)
3481  intvec* ivMat = MivMatrixOrder(iva);
3482
3483  int  mtmp, m=(*iva)[0];
3484
3485  for(i=ivMat->length(); i>=0; i--)
3486  {
3487    mtmp = (*ivMat)[i];
3488    if(mtmp <0)
3489    {
3490      mtmp = -mtmp;
3491    }
3492    if(mtmp > m)
3493    {
3494      m = mtmp;
3495    }
3496  }
3497
3498  // define the maximal total degree of polynomials of G
3499  mpz_t ndeg;
3500  mpz_init(ndeg);
3501
3502 // 12 Juli 03
3503#ifndef UPPER_BOUND
3504  mpz_set_si(ndeg, Trandegreebound(G)+1);
3505#else
3506  mpz_t ztmp;
3507  mpz_init(ztmp);
3508
3509  mpz_t maxdeg;
3510  mpz_init_set_si(maxdeg, Trandegreebound(G));
3511
3512  //ndeg = (2*maxdeg*maxdeg + (nV+1)*maxdeg)*m;//Kalkbrenner (1999)
3513  mpz_pow_ui(ztmp, maxdeg, 2);
3514  mpz_mul_ui(ztmp, ztmp, 2);
3515  mpz_mul_ui(maxdeg, maxdeg, nV+1);
3516  mpz_add(ndeg, ztmp, maxdeg);
3517  mpz_mul_ui(ndeg, ndeg, m);
3518
3519  mpz_clear(ztmp);
3520
3521  //PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m ");
3522  //Print("\n//         where d = %d, n = %d and bound = %d", maxdeg, nV, ndeg);
3523#endif //UPPER_BOUND
3524
3525#ifdef INVEPS_SMALL_IN_TRAN
3526  if(mpz_cmp_ui(ndeg, nV)>0 && nV > 3)
3527  {
3528    mpz_cdiv_q_ui(ndeg, ndeg, nV);
3529  }
3530 //PrintS("\n// choose the \"small\" inverse epsilon:");
3531 //mpz_out_str(stdout, 10, ndeg);
3532#endif
3533  mpz_t deg_tmp;
3534  mpz_init_set(deg_tmp, ndeg);
3535
3536  mpz_t *ivres=( mpz_t *) omAlloc(nV*sizeof(mpz_t));
3537  mpz_init_set_si(ivres[nV-1],1);
3538
3539  for(i=nV-2; i>=0; i--)
3540  {
3541    mpz_init_set(ivres[i], deg_tmp);
3542    mpz_mul(deg_tmp, deg_tmp, ndeg);
3543  }
3544
3545  mpz_t *ivtmp=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
3546  for(i=0; i<nV; i++)
3547  {
3548    mpz_init(ivtmp[i]);
3549  }
3550  mpz_t sing_int;
3551  mpz_init_set_ui(sing_int,  2147483647);
3552
3553  intvec* repr_vector = new intvec(nV);
3554
3555  // define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n
3556  for(i=0; i<nV; i++)
3557  {
3558    for(j=0; j<nV; j++)
3559    {
3560      if( (*ivMat)[i*nV+j] >= 0 )
3561      {
3562        mpz_mul_ui(ivres[i], ivres[i], (*ivMat)[i*nV+j]);
3563      }
3564      else
3565      {
3566        mpz_mul_ui(ivres[i], ivres[i], -(*ivMat)[i*nV+j]);
3567        mpz_neg(ivres[i], ivres[i]);
3568      }
3569      mpz_add(ivtmp[j], ivtmp[j], ivres[i]);
3570    }
3571  }
3572  delete ivMat;
3573
3574  int ntrue=0;
3575  for(i=0; i<nV; i++)
3576  {
3577    (*repr_vector)[i] = mpz_get_si(ivtmp[i]);
3578    if(mpz_cmp(ivtmp[i], sing_int)>=0)
3579    {
3580      ntrue++;
3581      if(Overflow_Error == FALSE)
3582      {
3583        Overflow_Error = TRUE;
3584
3585        PrintS("\n// ** OVERFLOW in \"Repr.Vector\": ");
3586        mpz_out_str( stdout, 10, ivtmp[i]);
3587        PrintS(" is greater than 2147483647 (max. integer representation)");
3588        Print("\n//  So vector[%d] := %d is wrong!!\n",i+1,(*repr_vector)[i]);
3589      }
3590    }
3591  }
3592  if(Overflow_Error == TRUE)
3593  {
3594    ivString(repr_vector, "repvector");
3595    Print("\n// %d element(s) of it are overflow!!", ntrue);
3596  }
3597
3598  if(Overflow_Error == FALSE)
3599    Overflow_Error=nError;
3600
3601  omFree(ivres);
3602  omFree(ivtmp);
3603
3604  mpz_clear(sing_int);
3605  mpz_clear(deg_tmp);
3606  mpz_clear(ndeg);
3607
3608  return repr_vector;
3609}
3610#endif
3611
3612//unused
3613#if 0
3614static intvec* TranPertVector_lp(ideal G)
3615{
3616  BOOLEAN nError = Overflow_Error;
3617  Overflow_Error = FALSE;
3618  // int j, nG = IDELEMS(G);
3619  int i;
3620  int nV = currRing->N;
3621
3622  // define the maximal total degree of polynomials of G
3623  mpz_t ndeg;
3624  mpz_init(ndeg);
3625
3626 // 12 Juli 03
3627#ifndef UPPER_BOUND
3628  mpz_set_si(ndeg, Trandegreebound(G)+1);
3629#else
3630  mpz_t ztmp;
3631  mpz_init(ztmp);
3632
3633  mpz_t maxdeg;
3634  mpz_init_set_si(maxdeg, Trandegreebound(G));
3635
3636  //ndeg = (2*maxdeg*maxdeg + (nV+1)*maxdeg);//Kalkbrenner (1999)
3637  mpz_pow_ui(ztmp, maxdeg, 2);
3638  mpz_mul_ui(ztmp, ztmp, 2);
3639  mpz_mul_ui(maxdeg, maxdeg, nV+1);
3640  mpz_add(ndeg, ztmp, maxdeg);
3641  // PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m ");
3642  // Print("\n//         where d = %d, n = %d and bound = %d",
3643  // mpz_get_si(maxdeg), nV, mpz_get_si(ndeg));
3644
3645 mpz_clear(ztmp);
3646
3647#endif
3648
3649#ifdef INVEPS_SMALL_IN_TRAN
3650 if(mpz_cmp_ui(ndeg, nV)>0 && nV > 3)
3651    mpz_cdiv_q_ui(ndeg, ndeg, nV);
3652
3653 //PrintS("\n// choose the \"small\" inverse epsilon:");
3654 // mpz_out_str(stdout, 10, ndeg);
3655#endif
3656
3657  mpz_t deg_tmp;
3658  mpz_init_set(deg_tmp, ndeg);
3659
3660  mpz_t *ivres=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
3661  mpz_init_set_si(ivres[nV-1], 1);
3662
3663  for(i=nV-2; i>=0; i--)
3664  {
3665    mpz_init_set(ivres[i], deg_tmp);
3666    mpz_mul(deg_tmp, deg_tmp, ndeg);
3667  }
3668
3669  mpz_t sing_int;
3670  mpz_init_set_ui(sing_int,  2147483647);
3671
3672  intvec* repr_vector = new intvec(nV);
3673  int ntrue=0;
3674  for(i=0; i<nV; i++)
3675  {
3676    (*repr_vector)[i] = mpz_get_si(ivres[i]);
3677
3678    if(mpz_cmp(ivres[i], sing_int)>=0)
3679    {
3680      ntrue++;
3681      if(Overflow_Error == FALSE)
3682      {
3683        Overflow_Error = TRUE;
3684        PrintS("\n// ** OVERFLOW in \"Repr.Vector\": ");
3685        mpz_out_str( stdout, 10, ivres[i]);
3686        PrintS(" is greater than 2147483647 (max. integer representation)");
3687        Print("\n//  So vector[%d] := %d is wrong!!\n",i+1,(*repr_vector)[i]);
3688      }
3689    }
3690  }
3691  if(Overflow_Error == TRUE)
3692  {
3693    ivString(repr_vector, "repvector");
3694    Print("\n// %d element(s) of it are overflow!!", ntrue);
3695  }
3696  if(Overflow_Error == FALSE)
3697    Overflow_Error = nError;
3698
3699  omFree(ivres);
3700
3701  mpz_clear(ndeg);
3702  mpz_clear(sing_int);
3703
3704  return repr_vector;
3705}
3706#endif
3707
3708//unused
3709#if 0
3710static intvec* RepresentationMatrix_Dp(ideal G, intvec* M)
3711{
3712  BOOLEAN nError = Overflow_Error;
3713  Overflow_Error = FALSE;
3714
3715  int i, j;
3716  int nV = currRing->N;
3717
3718  intvec* ivUnit = Mivdp(nV);
3719  int degtmp, maxdeg = 0;
3720
3721  for(i=IDELEMS(G)-1; i>=0; i--)
3722  {
3723    // find the maximal total degree of the polynomial G[i]
3724    degtmp = MwalkWeightDegree(G->m[i], ivUnit);
3725    if(degtmp > maxdeg)
3726      maxdeg = degtmp;
3727  }
3728
3729  mpz_t ztmp;
3730  mpz_init_set_si(ztmp, maxdeg);
3731  mpz_t *ivres=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
3732  mpz_init_set_si(ivres[nV-1], 1); // (*ivres)[nV-1] = 1;
3733
3734  for(i=nV-2; i>=0; i--)
3735  {
3736    mpz_init_set(ivres[i], ztmp); //(*ivres)[i] = ztmp;
3737    mpz_mul_ui(ztmp, ztmp, maxdeg); //ztmp *=maxdeg;
3738  }
3739
3740  mpz_t *ivtmp=(mpz_t*)omAlloc(nV*sizeof(mpz_t));
3741  for(i=0; i<nV; i++)
3742    mpz_init(ivtmp[i]);
3743
3744  // define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n
3745  for(i=0; i<nV; i++)
3746    for(j=0; j<nV; j++)
3747    {
3748      if((*M)[i*nV+j] < 0)
3749      {
3750        mpz_mul_ui(ztmp, ivres[i], -(*M)[i*nV+j]);
3751        mpz_neg(ztmp, ztmp);
3752      }
3753      else
3754        mpz_mul_ui(ztmp, ivres[i], (*M)[i*nV+j]);
3755
3756      mpz_add(ivtmp[j], ivtmp[j], ztmp);
3757    }
3758  delete ivres;
3759  mpz_t sing_int;
3760  mpz_init_set_ui(sing_int,  2147483647);
3761
3762  int ntrue=0;
3763  intvec* repvector = new intvec(nV);
3764  for(i=0; i<nV; i++)
3765  {
3766    (*repvector)[i] = mpz_get_si(ivtmp[i]);
3767    if(mpz_cmp(ivtmp[i], sing_int)>0)
3768    {
3769      ntrue++;
3770      if(Overflow_Error == FALSE)
3771      {
3772        Overflow_Error = TRUE;
3773        PrintS("\n// ** OVERFLOW in \"Repr.Matrix\": ");
3774        mpz_out_str( stdout, 10, ivtmp[i]);
3775        PrintS(" is greater than 2147483647 (max. integer representation)");
3776        Print("\n//  So vector[%d] := %d is wrong!!\n",i+1,(*repvector)[i]);
3777      }
3778    }
3779  }
3780  if(Overflow_Error == TRUE)
3781  {
3782    ivString(repvector, "repvector");
3783    Print("\n// %d element(s) of it are overflow!!", ntrue);
3784  }
3785
3786  if(Overflow_Error == FALSE)
3787    Overflow_Error = nError;
3788
3789  mpz_clear(sing_int);
3790  mpz_clear(ztmp);
3791  omFree(ivtmp);
3792  omFree(ivres);
3793  return repvector;
3794}
3795#endif
3796
3797/*****************************************************************************
3798 * The following subroutine is the implementation of our first improved      *
3799 * Groebner walk algorithm, i.e. the first altervative algorithm.            *
3800 * First we use the Grobner walk algorithm and then we call the changed      *
3801 * perturbation walk algorithm with decreased degree, if an intermediate     *
3802 * weight vector is equal to the current target weight vector.               *
3803 * This call will be only repeated until we get the wanted reduced Groebner  *
3804 * basis or n times, where n is the numbers of variables.                    *
3805 *****************************************************************************/
3806
3807//unused
3808#if 0
3809static int testnegintvec(intvec* v)
3810{
3811  int n = v->length();
3812  int i;
3813  for(i=0; i<n; i++)
3814  {
3815    if((*v)[i]<0)
3816    {
3817      return(1);
3818    }
3819  }
3820  return(0);
3821}
3822#endif
3823
3824// npwinc = 0, if curr_weight doesn't stay in the correct Groebner cone
3825static ideal Rec_LastGB(ideal G, intvec* curr_weight,
3826                        intvec* orig_target_weight, int tp_deg, int npwinc)
3827{
3828  BOOLEAN nError = Overflow_Error;
3829  Overflow_Error = FALSE;
3830  // BOOLEAN nOverflow_Error = FALSE;
3831
3832  clock_t tproc=0;
3833  clock_t tinput = clock();
3834
3835  int i,  nV = currRing->N;
3836  int nwalk=0, endwalks=0, nnwinC=1;
3837  int nlast = 0;
3838  ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG;
3839  ring newRing, oldRing, TargetRing;
3840  intvec* iv_M_lp;
3841  intvec* target_weight;
3842  intvec* ivNull = new intvec(nV); //define (0,...,0)
3843  ring EXXRing = currRing;
3844  //int NEG=0; //19 juni 03
3845  intvec* next_weight;
3846#ifndef  BUCHBERGER_ALG
3847  //08 Juli 03
3848  intvec* hilb_func;
3849#endif
3850  // to avoid (1,0,...,0) as the target vector
3851  intvec* last_omega = new intvec(nV);
3852  for(i=nV-1; i>0; i--)
3853    (*last_omega)[i] = 1;
3854  (*last_omega)[0] = 10000;
3855
3856  BOOLEAN isGB = FALSE;
3857
3858  // compute a pertubed weight vector of the target weight vector
3859  if(tp_deg > 1 && tp_deg <= nV)
3860  {
3861    ideal H0 = idHeadCC(G);
3862
3863    if (rParameter (currRing) != NULL)
3864    {
3865      DefRingParlp();
3866    }
3867    else
3868    {
3869      VMrDefaultlp();
3870    }
3871    TargetRing = currRing;
3872    ssG = idrMoveR(G,EXXRing,currRing);
3873
3874    ideal H0_tmp = idrMoveR(H0,EXXRing,currRing);
3875    ideal H1 = idHeadCC(ssG);
3876
3877    // Apply Lemma 2.2 in Collart et. al (1997) to check whether cone(k-1) is equal to cone(k)
3878    if(test_G_GB_walk(H0_tmp,H1)==1)
3879    {
3880      idDelete(&H0_tmp);
3881      idDelete(&H1);
3882      G = ssG;
3883      ssG = NULL;
3884      newRing = currRing;
3885      delete ivNull;
3886
3887      if(npwinc != 0)
3888      {
3889        goto LastGB_Finish;
3890      }
3891      else
3892      {
3893        isGB = TRUE;
3894        goto KSTD_Finish;
3895      }
3896    }
3897    idDelete(&H0_tmp);
3898    idDelete(&H1);
3899
3900    iv_M_lp = MivMatrixOrderlp(nV);
3901    target_weight  = MPertVectors(ssG, iv_M_lp, tp_deg);
3902    delete iv_M_lp;
3903    //PrintS("\n// Input is not GB!!");
3904    rChangeCurrRing(EXXRing);
3905    G = idrMoveR(ssG, TargetRing,currRing);
3906
3907    if(Overflow_Error == TRUE)
3908    {
3909      //nOverflow_Error = Overflow_Error;
3910      //NEG = 1;
3911      newRing = currRing;
3912      goto JUNI_STD;
3913    }
3914  }
3915
3916  while(1)
3917  {
3918    nwalk ++;
3919    nstep++;
3920
3921    if(nwalk==1)
3922    {
3923      goto FIRST_STEP;
3924    }
3925    to=clock();
3926    // compute an initial form ideal of <G> w.r.t. "curr_vector"
3927    Gomega = MwalkInitialForm(G, curr_weight);
3928    xtif=xtif+clock()-to;
3929
3930#ifndef  BUCHBERGER_ALG
3931    if(isNolVector(curr_weight) == 0)
3932    {
3933      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
3934    }
3935    else
3936    {
3937      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
3938    }
3939#endif // BUCHBERGER_ALG
3940
3941    oldRing = currRing;
3942
3943    // defiNe a new ring that its ordering is "(a(curr_weight),lp)
3944    if (rParameter(currRing) != NULL)
3945    {
3946      DefRingPar(curr_weight);
3947    }
3948    else
3949    {
3950      rChangeCurrRing(VMrDefault(curr_weight));
3951    }
3952    newRing = currRing;
3953    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
3954    to=clock();
3955    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
3956#ifdef  BUCHBERGER_ALG
3957    M = MstdhomCC(Gomega1);
3958#else
3959    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
3960    delete hilb_func;
3961#endif // BUCHBERGER_ALG
3962    xtstd=xtstd+clock()-to;
3963    // change the ring to oldRing
3964    rChangeCurrRing(oldRing);
3965    M1 =  idrMoveR(M, newRing,currRing);
3966    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
3967
3968     to=clock();
3969    // compute a reduced Groebner basis of <G> w.r.t. "newRing" by the lifting process
3970    F = MLifttwoIdeal(Gomega2, M1, G);
3971    xtlift=xtlift+clock()-to;
3972    idDelete(&M1);
3973    idDelete(&Gomega2);
3974    idDelete(&G);
3975
3976    // change the ring to newRing
3977    rChangeCurrRing(newRing);
3978    F1 = idrMoveR(F, oldRing,currRing);
3979
3980    to=clock();
3981    // reduce the Groebner basis <G> w.r.t. new ring
3982    G = kInterRedCC(F1, NULL);
3983    xtred=xtred+clock()-to;
3984    idDelete(&F1);
3985
3986    if(endwalks == 1)
3987    {
3988      break;
3989    }
3990  FIRST_STEP:
3991    to=clock();
3992    Overflow_Error = FALSE;
3993    // compute a next weight vector
3994    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
3995    xtnw=xtnw+clock()-to;
3996#ifdef PRINT_VECTORS
3997    MivString(curr_weight, target_weight, next_weight);
3998#endif
3999
4000    if(Overflow_Error == TRUE)
4001    {
4002      //PrintS("\n// ** The next vector does NOT stay in Cone!!\n");
4003#ifdef TEST_OVERFLOW
4004      goto  LastGB_Finish;
4005#endif
4006
4007      nnwinC = 0;
4008      if(tp_deg == nV)
4009      {
4010        nlast = 1;
4011      }
4012      delete next_weight;
4013      break;
4014    }
4015
4016    if(MivComp(next_weight, ivNull) == 1)
4017    {
4018      //newRing = currRing;
4019      delete next_weight;
4020      break;
4021    }
4022
4023    if(MivComp(next_weight, target_weight) == 1)
4024    {
4025      if(tp_deg == nV)
4026      {
4027        endwalks = 1;
4028      }
4029      else
4030      {
4031        // REC_LAST_GB_ALT2:
4032        //nOverflow_Error = Overflow_Error;
4033        tproc=tproc+clock()-tinput;
4034        /*
4035          Print("\n// takes %d steps and calls \"Rec_LastGB\" (%d):",
4036          nwalk, tp_deg+1);
4037        */
4038        G = Rec_LastGB(G,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
4039        newRing = currRing;
4040        delete next_weight;
4041        break;
4042      }
4043    }
4044
4045    for(i=nV-1; i>=0; i--)
4046    {
4047      (*curr_weight)[i] = (*next_weight)[i];
4048    }
4049    delete next_weight;
4050  }//while
4051
4052  delete ivNull;
4053
4054  if(tp_deg != nV)
4055  {
4056    newRing = currRing;
4057
4058    if (rParameter(currRing) != NULL)
4059    {
4060      DefRingParlp();
4061    }
4062    else
4063    {
4064      VMrDefaultlp();
4065    }
4066    F1 = idrMoveR(G, newRing,currRing);
4067
4068    if(nnwinC == 0 || test_w_in_ConeCC(F1, target_weight) != 1 )
4069    {
4070      // nOverflow_Error = Overflow_Error;
4071      //Print("\n//  takes %d steps and calls \"Rec_LastGB (%d):", tp_deg+1);
4072      tproc=tproc+clock()-tinput;
4073      F1 = Rec_LastGB(F1,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
4074    }
4075    delete target_weight;
4076
4077    TargetRing = currRing;
4078    rChangeCurrRing(EXXRing);
4079    result = idrMoveR(F1, TargetRing,currRing);
4080  }
4081  else
4082  {
4083    if(nlast == 1)
4084    {
4085      JUNI_STD:
4086
4087      newRing = currRing;
4088      if (rParameter(currRing) != NULL)
4089      {
4090        DefRingParlp();
4091      }
4092      else
4093      {
4094        VMrDefaultlp();
4095      }
4096      KSTD_Finish:
4097      if(isGB == FALSE)
4098      {
4099        F1 = idrMoveR(G, newRing,currRing);
4100      }
4101      else
4102      {
4103        F1 = G;
4104      }
4105      to=clock();
4106      // Print("\n// apply the Buchberger's alg in ring = %s",rString(currRing));
4107      // idElements(F1, "F1");
4108      G = MstdCC(F1);
4109      xtextra=xtextra+clock()-to;
4110
4111
4112      idDelete(&F1);
4113      newRing = currRing;
4114    }
4115
4116    LastGB_Finish:
4117    rChangeCurrRing(EXXRing);
4118    result = idrMoveR(G, newRing,currRing);
4119  }
4120
4121  if(Overflow_Error == FALSE)
4122    {
4123    Overflow_Error=nError;
4124    }
4125// Print("\n// \"Rec_LastGB\" (%d) took %d steps and %.2f sec.Overflow_Error (%d)", tp_deg, nwalk, ((double) tproc)/1000000, nOverflow_Error);
4126  return(result);
4127}
4128
4129/* The following subroutine is the implementation of our second improved
4130   Groebner walk algorithm, i.e. the second altervative algorithm.
4131   First we use the Grobner walk algorithm and then we call the changed
4132   perturbation walk algorithm with increased degree, if an intermediate
4133   weight vector is equal to the current target weight vector.
4134   This call will be only repeated until we get the wanted reduced Groebner
4135   basis or n times, where n is the numbers of variables.
4136*/
4137
4138/******************************
4139 * walk + recursive LastGB    *
4140 ******************************/
4141ideal MAltwalk2(ideal Go, intvec* curr_weight, intvec* target_weight)
4142{
4143  Set_Error(FALSE);
4144  Overflow_Error = FALSE;
4145  //BOOLEAN nOverflow_Error = FALSE;
4146  //Print("// pSetm_Error = (%d)", ErrorCheck());
4147
4148  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0;
4149  xftinput = clock();
4150  clock_t tostd, tproc;
4151
4152  nstep = 0;
4153  int i, nV = currRing->N;
4154  int nwalk=0, endwalks=0;
4155  // int nhilb = 1;
4156  ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G;
4157  //ideal  G1;
4158  //ring endRing;
4159  ring newRing, oldRing;
4160  intvec* ivNull = new intvec(nV);
4161  intvec* next_weight;
4162#if 0
4163  intvec* extra_curr_weight = new intvec(nV);
4164#endif
4165  //intvec* hilb_func;
4166  intvec* exivlp = Mivlp(nV);
4167
4168  ring XXRing = currRing;
4169
4170  //Print("\n// ring r_input = %s;", rString(currRing));
4171  to = clock();
4172  /* compute the reduced Groebner basis of the given ideal w.r.t.
4173     a "fast" monomial order, e.g. degree reverse lex. order (dp) */
4174  G = MstdCC(Go);
4175  tostd=clock()-to;
4176
4177  /*
4178  Print("\n// Computation of the first std took = %.2f sec",
4179        ((double) tostd)/1000000);
4180  */
4181  if(currRing->order[0] == ringorder_a)
4182  {
4183    goto NEXT_VECTOR;
4184  }
4185  while(1)
4186  {
4187    nwalk ++;
4188    nstep ++;
4189    to = clock();
4190    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
4191    Gomega = MwalkInitialForm(G, curr_weight);
4192    xtif=xtif+clock()-to;
4193#if 0
4194    if(Overflow_Error == TRUE)
4195    {
4196      for(i=nV-1; i>=0; i--)
4197        (*curr_weight)[i] = (*extra_curr_weight)[i];
4198      delete extra_curr_weight;
4199      goto LAST_GB_ALT2;
4200    }
4201#endif
4202    oldRing = currRing;
4203
4204    /* define a new ring that its ordering is "(a(curr_weight),lp) */
4205    if (rParameter(currRing) != NULL)
4206    {
4207      DefRingPar(curr_weight);
4208    }
4209    else
4210    {
4211      rChangeCurrRing(VMrDefault(curr_weight));
4212    }
4213    newRing = currRing;
4214    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
4215    to = clock();
4216    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
4217    M = MstdhomCC(Gomega1);
4218    xtstd=xtstd+clock()-to;
4219    /* change the ring to oldRing */
4220    rChangeCurrRing(oldRing);
4221    M1 =  idrMoveR(M, newRing,currRing);
4222    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
4223
4224    to = clock();
4225    /* compute the reduced Groebner basis of <G> w.r.t. "newRing"
4226       by the liftig process */
4227    F = MLifttwoIdeal(Gomega2, M1, G);
4228    xtlift=xtlift+clock()-to;
4229    idDelete(&M1);
4230    idDelete(&Gomega2);
4231    idDelete(&G);
4232
4233    /* change the ring to newRing */
4234    rChangeCurrRing(newRing);
4235    F1 = idrMoveR(F, oldRing,currRing);
4236
4237    to = clock();
4238    /* reduce the Groebner basis <G> w.r.t. newRing */
4239    G = kInterRedCC(F1, NULL);
4240    xtred=xtred+clock()-to;
4241    idDelete(&F1);
4242
4243    if(endwalks == 1)
4244      break;
4245
4246  NEXT_VECTOR:
4247    to = clock();
4248    /* compute a next weight vector */
4249    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
4250    xtnw=xtnw+clock()-to;
4251#ifdef PRINT_VECTORS
4252    MivString(curr_weight, target_weight, next_weight);
4253#endif
4254
4255    if(Overflow_Error == TRUE)
4256    {
4257      /*
4258        ivString(next_weight, "omega");
4259        PrintS("\n// ** The weight vector does NOT stay in Cone!!\n");
4260      */
4261#ifdef TEST_OVERFLOW
4262      goto  TEST_OVERFLOW_OI;
4263#endif
4264
4265      newRing = currRing;
4266      if (rParameter(currRing) != NULL)
4267      {
4268        DefRingPar(target_weight);
4269      }
4270      else
4271      {
4272        rChangeCurrRing(VMrDefault(target_weight)); // Aenderung
4273      }
4274      F1 = idrMoveR(G, newRing,currRing);
4275      G = MstdCC(F1);
4276      idDelete(&F1);
4277      newRing = currRing;
4278      break;
4279    }
4280
4281    if(MivComp(next_weight, ivNull) == 1)
4282    {
4283      newRing = currRing;
4284      delete next_weight;
4285      break;
4286    }
4287
4288    if(MivComp(next_weight, target_weight) == 1)
4289    {
4290      if(MivSame(target_weight, exivlp)==1)
4291      {
4292     // LAST_GB_ALT2:
4293        //nOverflow_Error = Overflow_Error;
4294        tproc = clock()-xftinput;
4295        //Print("\n// takes %d steps and calls the recursion of level 2:",  nwalk);
4296        /* call the changed perturbation walk algorithm with degree 2 */
4297        G = Rec_LastGB(G, curr_weight, target_weight, 2,1);
4298        newRing = currRing;
4299        delete next_weight;
4300        break;
4301      }
4302      endwalks = 1;
4303    }
4304
4305    for(i=nV-1; i>=0; i--)
4306    {
4307      //(*extra_curr_weight)[i] = (*curr_weight)[i];
4308      (*curr_weight)[i] = (*next_weight)[i];
4309    }
4310    delete next_weight;
4311  }
4312#ifdef TEST_OVERFLOW
4313 TEST_OVERFLOW_OI:
4314#endif
4315  rChangeCurrRing(XXRing);
4316  G = idrMoveR(G, newRing,currRing);
4317  delete ivNull;
4318  delete exivlp;
4319
4320#ifdef TIME_TEST
4321 // Print("\n// \"Main procedure\"  took %d steps dnd %.2f sec. Overflow_Error (%d)", nwalk, ((double) tproc)/1000000, nOverflow_Error);
4322
4323  TimeStringFractal(xftinput, tostd, xtif, xtstd, xtextra,xtlift, xtred,xtnw);
4324
4325  Print("\n// pSetm_Error = (%d)", ErrorCheck());
4326  //Print("\n// Overflow_Error? (%d)", nOverflow_Error);
4327  Print("\n// Awalk2 took %d steps!!", nstep);
4328#endif
4329
4330  return(G);
4331}
4332
4333
4334/**************************************
4335 * perturb the matrix order of  "lex" *
4336 **************************************/
4337static intvec* NewVectorlp(ideal I)
4338{
4339  int nV = currRing->N;
4340  intvec* iv_wlp =  MivMatrixOrderlp(nV);
4341  intvec* result = Mfpertvector(I, iv_wlp);
4342  delete iv_wlp;
4343  return result;
4344}
4345
4346int ngleich;
4347intvec* Xsigma;
4348intvec* Xtau;
4349int xn;
4350intvec* Xivinput;
4351intvec* Xivlp;
4352
4353
4354/********************************
4355 * compute a next weight vector *
4356 ********************************/
4357static intvec* MWalkRandomNextWeight(ideal G, intvec* curr_weight,
4358               intvec* target_weight, int weight_rad, int pert_deg)
4359{
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);
4371  intvec* next_weight22 = new intvec(nV);
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)
4415      {
4416        for(i=0; i<nV; i++)
4417        {
4418          (*next_weight22)[i] = (*next_weight2)[i];
4419        }
4420        i = 0;
4421        // reduce the size of the maximal entry of the vector
4422        while(test_w_in_ConeCC(G,next_weight22) == 1)
4423        {
4424          (*next_weight2)[i] = (*next_weight22)[i];
4425          i = MivAbsMaxArg(next_weight22);
4426          (*next_weight22)[i] = floor(0.1*(*next_weight22)[i] + 0.5);
4427        }
4428        delete next_weight22;
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  {
4443    ideal G_test1 = MwalkInitialForm(G, next_weight1);
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))
4447      {
4448        for(i=0; i<nV; i++)
4449        {
4450          (*result)[i] = (*next_weight2)[i];
4451        }
4452      }
4453      else
4454      {
4455        for(i=0; i<nV; i++)
4456        {
4457          (*result)[i] = (*next_weight1)[i];
4458        }
4459      }   
4460    }
4461    else
4462    {
4463      if(G_test2->m[0] != NULL && maxlengthpoly(G_test2) < maxlengthpoly(G_test))
4464      {
4465        for(i=0; i<nV; i++)
4466        {
4467          (*result)[i] = (*next_weight2)[i];
4468        }
4469      }
4470      else
4471      {
4472        for(i=0; i<nV; i++)
4473        {
4474          (*result)[i] = (*next_weight)[i];
4475        }
4476      }
4477    }
4478    idDelete(&G_test1);
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      }
4489    }
4490    else
4491    {
4492      for(i=0; i<nV; i++)
4493      {
4494        (*result)[i] = (*next_weight)[i];
4495      }
4496    }
4497  }
4498
4499  idDelete(&G_test);
4500  idDelete(&G_test2);
4501  if(test_w_in_ConeCC(G, result) == 1)
4502  {
4503    delete next_weight2;
4504    delete next_weight;
4505    delete next_weight1;
4506    return result;
4507  }
4508  else
4509  {
4510    delete result;
4511    delete next_weight2;
4512    delete next_weight1;
4513    return next_weight;
4514  }
4515}
4516
4517
4518/***************************************************************************
4519 * The procedur REC_GB_Mwalk computes a GB for <G> w.r.t. the weight order *
4520 * otw, where G is a reduced GB w.r.t. the weight order cw.                *
4521 * The new procedur Mwalk calls REC_GB.                                    *
4522 ***************************************************************************/
4523static ideal REC_GB_Mwalk(ideal G, intvec* curr_weight, intvec* orig_target_weight,
4524                          int tp_deg, int npwinc)
4525{
4526  BOOLEAN nError = Overflow_Error;
4527  Overflow_Error = FALSE;
4528
4529  int i,  nV = currRing->N;
4530  int nwalk=0, endwalks=0, nnwinC=1, nlast = 0;
4531  ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG;
4532  ring newRing, oldRing, TargetRing;
4533  intvec* target_weight;
4534  intvec* ivNull = new intvec(nV);
4535#ifndef BUCHBERGER_ALG
4536  intvec* hilb_func;
4537  // to avoid (1,0,...,0) as the target vector
4538  intvec* last_omega = new intvec(nV);
4539  for(i=nV-1; i>0; i--)
4540  {
4541    (*last_omega)[i] = 1;
4542  }
4543  (*last_omega)[0] = 10000;
4544#endif
4545  BOOLEAN isGB = FALSE;
4546
4547  ring EXXRing = currRing;
4548
4549  // compute a pertubed weight vector of the target weight vector
4550  if(tp_deg > 1 && tp_deg <= nV)
4551  {
4552    ideal H0 = idHeadCC(G);
4553    if (rParameter(currRing) != NULL)
4554    {
4555      DefRingPar(orig_target_weight);
4556    }
4557    else
4558    {
4559      rChangeCurrRing(VMrDefault(orig_target_weight));
4560    }
4561    TargetRing = currRing;
4562    ssG = idrMoveR(G,EXXRing,currRing);
4563
4564    ideal H0_tmp = idrMoveR(H0,EXXRing,currRing);
4565    ideal H1 = idHeadCC(ssG);
4566    id_Delete(&H0,EXXRing);
4567
4568    if(test_G_GB_walk(H0_tmp,H1)==1)
4569    {
4570      //Print("\n//REC_GB_Mwalk: input in %d-th recursive is a GB!\n",tp_deg);
4571      idDelete(&H0_tmp);
4572      idDelete(&H1);
4573      G = ssG;
4574      ssG = NULL;
4575      newRing = currRing;
4576      delete ivNull;
4577      if(npwinc == 0)
4578      {
4579        isGB = TRUE;
4580        goto KSTD_Finish;
4581      }
4582      else
4583      {
4584        goto LastGB_Finish;
4585      }
4586    }
4587    idDelete(&H0_tmp);
4588    idDelete(&H1);
4589
4590    target_weight  = MPertVectors(ssG, MivMatrixOrder(orig_target_weight), tp_deg);
4591
4592    rChangeCurrRing(EXXRing);
4593    G = idrMoveR(ssG, TargetRing,currRing);
4594  }
4595
4596  while(1)
4597  {
4598    nwalk ++;
4599    nstep++;
4600    if(nwalk == 1)
4601    {
4602      goto NEXT_STEP;
4603    }
4604    //Print("\n//REC_GB_Mwalk: Entering the %d-th step in the %d-th recursive:\n",nwalk,tp_deg);
4605    to = clock();
4606    // compute an initial form ideal of <G> w.r.t. "curr_vector"
4607    Gomega = MwalkInitialForm(G, curr_weight);
4608    xtif = xtif + clock()-to;
4609
4610#ifndef  BUCHBERGER_ALG
4611    if(isNolVector(curr_weight) == 0)
4612    {
4613      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
4614    }
4615    else
4616    {
4617      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
4618    }
4619#endif
4620
4621    oldRing = currRing;
4622
4623    // define a new ring with ordering "(a(curr_weight),lp)
4624    if (rParameter(currRing) != NULL)
4625    {
4626      DefRingPar(curr_weight);
4627    }
4628    else
4629    {
4630      rChangeCurrRing(VMrDefault(curr_weight));
4631    }
4632    newRing = currRing;
4633    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
4634
4635    to = clock();
4636    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
4637#ifdef  BUCHBERGER_ALG
4638    M = MstdhomCC(Gomega1);
4639#else
4640    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
4641    delete hilb_func;
4642#endif
4643    xtstd = xtstd + clock() - to;
4644
4645    // change the ring to oldRing
4646    rChangeCurrRing(oldRing);
4647
4648    M1 =  idrMoveR(M, newRing,currRing);
4649    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
4650
4651    to = clock();
4652    F = MLifttwoIdeal(Gomega2, M1, G);
4653    xtlift = xtlift + clock() -to;
4654
4655    idDelete(&M1);
4656    idDelete(&Gomega2);
4657    idDelete(&G);
4658
4659
4660    // change the ring to newRing
4661    rChangeCurrRing(newRing);
4662    F1 = idrMoveR(F, oldRing,currRing);
4663
4664    to = clock();
4665    // reduce the Groebner basis <G> w.r.t. new ring
4666    G = kInterRedCC(F1, NULL);
4667    xtred = xtred + clock() -to;
4668
4669    idDelete(&F1);
4670
4671    if(endwalks == 1)
4672    {
4673      break;
4674    }
4675  NEXT_STEP:
4676    to = clock();
4677    // compute a next weight vector
4678    intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
4679
4680
4681    xtnw = xtnw + clock() - to;
4682
4683#ifdef PRINT_VECTORS
4684    MivString(curr_weight, target_weight, next_weight);
4685#endif
4686
4687    if(Overflow_Error == TRUE)
4688    {
4689      //PrintS("\n//REC_GB_Mwalk: The computed vector does NOT stay in the correct cone!!\n");
4690      nnwinC = 0;
4691      if(tp_deg == nV)
4692      {
4693        nlast = 1;
4694      }
4695      delete next_weight;
4696      break;
4697    }
4698    if(MivComp(next_weight, ivNull) == 1)
4699    {
4700      newRing = currRing;
4701      delete next_weight;
4702      break;
4703    }
4704
4705    if(MivComp(next_weight, target_weight) == 1)
4706    {
4707      if(tp_deg == nV)
4708      {
4709        endwalks = 1;
4710      }
4711      else
4712      {
4713        G = REC_GB_Mwalk(G,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
4714        newRing = currRing;
4715        delete next_weight;
4716        break;
4717      }
4718    }
4719
4720    for(i=nV-1; i>=0; i--)
4721    {
4722      (*curr_weight)[i] = (*next_weight)[i];
4723    }
4724    delete next_weight;
4725  }
4726
4727  delete ivNull;
4728
4729  if(tp_deg != nV)
4730  {
4731    newRing = currRing;
4732
4733    if (rParameter(currRing) != NULL)
4734    {
4735      DefRingPar(orig_target_weight);
4736    }
4737    else
4738    {
4739      rChangeCurrRing(VMrDefault(orig_target_weight));
4740    }
4741    F1 = idrMoveR(G, newRing,currRing);
4742
4743    if(nnwinC == 0)
4744    {
4745      F1 = REC_GB_Mwalk(F1,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
4746    }
4747    else
4748    {
4749      if(test_w_in_ConeCC(F1, target_weight) != 1)
4750      {
4751        F1 = REC_GB_Mwalk(F1,curr_weight, orig_target_weight,tp_deg+1,nnwinC);
4752      }
4753    }
4754    delete target_weight;
4755
4756    TargetRing = currRing;
4757    rChangeCurrRing(EXXRing);
4758    result = idrMoveR(F1, TargetRing,currRing);
4759  }
4760  else
4761  {
4762    if(nlast == 1)
4763    {
4764      if (rParameter(currRing) != NULL)
4765      {
4766        DefRingPar(orig_target_weight);
4767      }
4768      else
4769      {
4770        rChangeCurrRing(VMrDefault(orig_target_weight));
4771      }
4772    KSTD_Finish:
4773      if(isGB == FALSE)
4774      {
4775        F1 = idrMoveR(G, newRing,currRing);
4776      }
4777      else
4778      {
4779        F1 = G;
4780      }
4781      to=clock();
4782      // apply Buchberger alg to compute a red. GB of F1
4783      G = MstdCC(F1);
4784      xtextra=clock()-to;
4785      idDelete(&F1);
4786      newRing = currRing;
4787    }
4788
4789  LastGB_Finish:
4790    rChangeCurrRing(EXXRing);
4791    result = idrMoveR(G, newRing,currRing);
4792  }
4793
4794  if(Overflow_Error == FALSE)
4795    {
4796    Overflow_Error = nError;
4797    }
4798#ifndef BUCHBERGER_ALG
4799  delete last_omega;
4800#endif
4801  return(result);
4802}
4803
4804
4805// THE NEW GROEBNER WALK ALGORITHM
4806// Groebnerwalk with a recursive "second" alternative GW, called REC_GB_Mwalk that only computes the last reduced GB
4807ideal MwalkAlt(ideal Go, intvec* curr_weight, intvec* target_weight)
4808{
4809  Set_Error(FALSE);
4810  Overflow_Error = FALSE;
4811  //Print("// pSetm_Error = (%d)", ErrorCheck());
4812
4813  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
4814  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
4815  tinput = clock();
4816  clock_t tim;
4817  nstep=0;
4818  int i;
4819  int nV = currRing->N;
4820  int nwalk=0;
4821  int endwalks=0;
4822
4823  ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G;
4824  //ideal G1;
4825  //ring endRing;
4826  ring newRing, oldRing;
4827  intvec* ivNull = new intvec(nV);
4828  intvec* exivlp = Mivlp(nV);
4829#ifndef BUCHBERGER_ALG
4830  intvec* hilb_func;
4831#endif
4832  intvec* tmp_weight = new intvec(nV);
4833  for(i=nV-1; i>=0; i--)
4834    (*tmp_weight)[i] = (*curr_weight)[i];
4835
4836   // to avoid (1,0,...,0) as the target vector
4837  intvec* last_omega = new intvec(nV);
4838  for(i=nV-1; i>0; i--)
4839    (*last_omega)[i] = 1;
4840  (*last_omega)[0] = 10000;
4841
4842  ring XXRing = currRing;
4843
4844  to = clock();
4845  // the monomial ordering of this current ring would be "dp"
4846  G = MstdCC(Go);
4847  tostd = clock()-to;
4848
4849  if(currRing->order[0] == ringorder_a)
4850    goto NEXT_VECTOR;
4851
4852  while(1)
4853  {
4854    nwalk ++;
4855    nstep ++;
4856    to = clock();
4857    // compute an initial form ideal of <G> w.r.t. "curr_vector"
4858    Gomega = MwalkInitialForm(G, curr_weight);
4859    tif = tif + clock()-to;
4860    oldRing = currRing;
4861
4862    if(endwalks == 1)
4863    {
4864      /* compute a reduced Groebner basis of Gomega w.r.t. >>_cw by
4865         the recursive changed perturbation walk alg. */
4866      tim = clock();
4867      /*
4868        Print("\n// **** Groebnerwalk took %d steps and ", nwalk);
4869        PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:");
4870        idElements(Gomega, "G_omega");
4871      */
4872
4873      if(MivSame(exivlp, target_weight)==1)
4874        M = REC_GB_Mwalk(idCopy(Gomega), tmp_weight, curr_weight, 2,1);
4875      else
4876        goto NORMAL_GW;
4877      /*
4878        Print("\n//  time for the last std(Gw)  = %.2f sec",
4879        ((double) (clock()-tim)/1000000));
4880        PrintS("\n// ***************************************************\n");
4881      */
4882#ifdef CHECK_IDEAL_MWALK
4883      idElements(Gomega, "G_omega");
4884      headidString(Gomega, "Gw");
4885      idElements(M, "M");
4886      //headidString(M, "M");
4887#endif
4888      to = clock();
4889      F = MLifttwoIdeal(Gomega, M, G);
4890      xtlift = xtlift + clock() - to;
4891
4892      idDelete(&Gomega);
4893      idDelete(&M);
4894      idDelete(&G);
4895
4896      oldRing = currRing;
4897
4898      // create a new ring newRing
4899       if (rParameter(currRing) != NULL)
4900       {
4901         DefRingPar(curr_weight);
4902       }
4903       else
4904       {
4905         rChangeCurrRing(VMrDefault(curr_weight));
4906       }
4907      newRing = currRing;
4908      F1 = idrMoveR(F, oldRing,currRing);
4909    }
4910    else
4911    {
4912    NORMAL_GW:
4913#ifndef  BUCHBERGER_ALG
4914      if(isNolVector(curr_weight) == 0)
4915      {
4916        hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
4917      }
4918      else
4919      {
4920        hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
4921      }
4922#endif // BUCHBERGER_ALG
4923
4924      // define a new ring that its ordering is "(a(curr_weight),lp)
4925      if (rParameter(currRing) != NULL)
4926      {
4927        DefRingPar(curr_weight);
4928      }
4929      else
4930      {
4931        rChangeCurrRing(VMrDefault(curr_weight));
4932      }
4933      newRing = currRing;
4934      Gomega1 = idrMoveR(Gomega, oldRing,currRing);
4935
4936      to = clock();
4937      // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
4938#ifdef  BUCHBERGER_ALG
4939      M = MstdhomCC(Gomega1);
4940#else
4941      M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
4942      delete hilb_func;
4943#endif
4944      tstd = tstd + clock() - to;
4945
4946      // change the ring to oldRing
4947      rChangeCurrRing(oldRing);
4948      M1 =  idrMoveR(M, newRing,currRing);
4949      Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
4950
4951      to = clock();
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.
4955      F = MLifttwoIdeal(Gomega2, M1, G);
4956      tlift = tlift + clock() - to;
4957
4958      idDelete(&M1);
4959      idDelete(&Gomega2);
4960      idDelete(&G);
4961
4962      // change the ring to newRing
4963      rChangeCurrRing(newRing);
4964      F1 = idrMoveR(F, oldRing,currRing);
4965    }
4966
4967    to = clock();
4968    // reduce the Groebner basis <G> w.r.t. new ring
4969    G = kInterRedCC(F1, NULL);
4970    if(endwalks != 1)
4971    {
4972      tred = tred + clock() - to;
4973    }
4974    else
4975    {
4976      xtred = xtred + clock() - to;
4977    }
4978    idDelete(&F1);
4979    if(endwalks == 1)
4980    {
4981      break;
4982    }
4983  NEXT_VECTOR:
4984    to = clock();
4985    // compute a next weight vector
4986    intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G);
4987    tnw = tnw + clock() - to;
4988#ifdef PRINT_VECTORS
4989    MivString(curr_weight, target_weight, next_weight);
4990#endif
4991
4992    //if(test_w_in_ConeCC(G, next_weight) != 1)
4993    if(Overflow_Error == TRUE)
4994    {
4995      newRing = currRing;
4996      PrintS("\n// ** The computed vector does NOT stay in Cone!!\n");
4997
4998      if (rParameter(currRing) != NULL)
4999      {
5000        DefRingPar(target_weight);
5001      }
5002      else
5003      {
5004        rChangeCurrRing(VMrDefault(target_weight));
5005      }
5006      F1 = idrMoveR(G, newRing,currRing);
5007      G = MstdCC(F1);
5008      idDelete(&F1);
5009
5010      newRing = currRing;
5011      break;
5012    }
5013
5014    if(MivComp(next_weight, ivNull) == 1)
5015    {
5016      newRing = currRing;
5017      delete next_weight;
5018      break;
5019    }
5020    if(MivComp(next_weight, target_weight) == 1)
5021    {
5022      endwalks = 1;
5023    }
5024    for(i=nV-1; i>=0; i--)
5025    {
5026      (*tmp_weight)[i] = (*curr_weight)[i];
5027      (*curr_weight)[i] = (*next_weight)[i];
5028    }
5029    delete next_weight;
5030  }
5031  rChangeCurrRing(XXRing);
5032  G = idrMoveR(G, newRing,currRing);
5033
5034  delete tmp_weight;
5035  delete ivNull;
5036  delete exivlp;
5037
5038#ifdef TIME_TEST
5039  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
5040
5041  Print("\n// pSetm_Error = (%d)", ErrorCheck());
5042  Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
5043#endif
5044  return(G);
5045}
5046
5047
5048/*******************************
5049 * THE GROEBNER WALK ALGORITHM *
5050 *******************************/
5051ideal Mwalk(ideal Go, intvec* orig_M, intvec* target_M,
5052            ring baseRing, int reduction, int printout)
5053{
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  }
5061  Set_Error(FALSE);
5062  Overflow_Error = FALSE;
5063  //BOOLEAN endwalks = FALSE;
5064#ifdef TIME_TEST
5065  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
5066  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
5067  tinput = clock();
5068  clock_t tim;
5069#endif
5070  nstep=0;
5071  int i,nwalk;
5072  int nV = baseRing->N;
5073
5074  ideal Gomega, M, F, FF, Gomega1, Gomega2, M1;
5075  ring newRing;
5076  ring XXRing = baseRing;
5077  ring targetRing;
5078  intvec* ivNull = new intvec(nV);
5079  intvec* curr_weight = new intvec(nV);
5080  intvec* target_weight = new intvec(nV);
5081  intvec* exivlp = Mivlp(nV);
5082/*
5083  intvec* tmp_weight = new intvec(nV);
5084  for(i=0; i<nV; i++)
5085  {
5086    (*tmp_weight)[i] = (*orig_M)[i];
5087  }
5088*/
5089  for(i=0; i<nV; i++)
5090  {
5091    (*curr_weight)[i] = (*orig_M)[i];
5092    (*target_weight)[i] = (*target_M)[i];
5093  }
5094#ifndef BUCHBERGER_ALG
5095  intvec* hilb_func;
5096   // to avoid (1,0,...,0) as the target vector
5097  intvec* last_omega = new intvec(nV);
5098  for(i=nV-1; i>0; i--)
5099  {
5100    (*last_omega)[i] = 1;
5101  }
5102  (*last_omega)[0] = 10000;
5103#endif
5104  rComplete(currRing);
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
5132#ifdef TIME_TEST
5133  to = clock();
5134#endif
5135
5136  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
5137
5138#ifdef TIME_TEST
5139  tostd = clock()-to;
5140#endif
5141
5142  baseRing = currRing;
5143  nwalk = 0;
5144
5145  while(1)
5146  {
5147    nwalk ++;
5148    nstep ++;
5149
5150#ifdef TIME_TEST
5151    to = clock();
5152#endif
5153    // compute an initial form ideal of <G> w.r.t. "curr_vector"
5154    Gomega = MwalkInitialForm(G, curr_weight);
5155#ifdef TIME_TEST
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
5178#ifndef  BUCHBERGER_ALG
5179    if(isNolVector(curr_weight) == 0)
5180    {
5181      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
5182    }   
5183    else
5184    {
5185      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
5186    }
5187#endif
5188
5189    if(nwalk == 1)
5190    {
5191      if(orig_M->length() == nV)
5192      {
5193        // define a new ring with ordering "(a(curr_weight),lp)
5194        newRing = VMrDefault(curr_weight);
5195      }
5196      else
5197      {
5198        newRing = VMatrDefault(orig_M);
5199      }
5200    }
5201    else
5202    {
5203     if(target_M->length() == nV)
5204     {
5205       //define a new ring with ordering "(a(curr_weight),lp)"
5206       newRing = VMrDefault(curr_weight);
5207     }
5208     else
5209     {
5210       //define a new ring with matrix ordering
5211       newRing = VMatrRefine(target_M,curr_weight);
5212     }
5213    }
5214    rChangeCurrRing(newRing);
5215    Gomega1 = idrMoveR(Gomega, baseRing,currRing);
5216    idDelete(&Gomega);
5217    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
5218#ifdef TIME_TEST
5219    to = clock();
5220#endif
5221#ifndef  BUCHBERGER_ALG
5222    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
5223    delete hilb_func;
5224#else
5225    M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
5226#endif
5227#ifdef TIME_TEST
5228    tstd = tstd + clock() - to;
5229#endif
5230    idSkipZeroes(M);
5231//#ifdef CHECK_IDEAL_MWALK
5232    if(printout > 2)
5233    {
5234      idString(M, "//** Mwalk: M");
5235    }
5236//#endif
5237    //change the ring to baseRing
5238    rChangeCurrRing(baseRing);
5239    M1 =  idrMoveR(M, newRing,currRing);
5240    idDelete(&M);
5241    Gomega2 = idrMoveR(Gomega1, newRing,currRing);
5242    idDelete(&Gomega1);
5243#ifdef TIME_TEST
5244    to = clock();
5245#endif
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
5248    F = MLifttwoIdeal(Gomega2, M1, G);
5249
5250#ifdef TIME_TEST
5251    tlift = tlift + clock() - to;
5252#endif
5253//#ifdef CHECK_IDEAL_MWALK
5254    if(printout > 2)
5255    {
5256      idString(F, "//** Mwalk: F");
5257    }
5258//#endif
5259    idDelete(&Gomega2);
5260    idDelete(&M1);
5261
5262    rChangeCurrRing(newRing); // change the ring to newRing
5263    G = idrMoveR(F,baseRing,currRing);
5264    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);
5286    baseRing = currRing;
5287
5288/*
5289#ifdef TIME_TEST
5290    to = clock();
5291#endif
5292
5293#ifdef TIME_TEST
5294    tstd = tstd + clock() - to;
5295#endif
5296*/
5297
5298
5299#ifdef TIME_TEST
5300    to = clock();
5301#endif
5302    NEXT_VECTOR:
5303    intvec* next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
5304#ifdef TIME_TEST
5305    tnw = tnw + clock() - to;
5306#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
5322      Gomega = MwalkInitialForm(G, curr_weight); // compute an initial form ideal of <G> w.r.t. "curr_vector"
5323      if(target_M->length() == nV)
5324      {
5325        newRing = VMrDefault(target_weight); // define a new ring with ordering "(a(curr_weight),lp)
5326      }
5327      else
5328      {
5329        newRing = VMatrDefault(target_M);
5330      }
5331      rChangeCurrRing(newRing);
5332      Gomega1 = idrMoveR(Gomega, baseRing,currRing);
5333      idDelete(&Gomega);
5334//#ifdef CHECK_IDEAL_MWALK
5335      if(printout > 1)
5336      {
5337        idString(Gomega1, "//** Mwalk: Gomega");
5338      }
5339      PrintS("\n //** Mwalk: kStd(Gomega)");
5340//#endif
5341      M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
5342//#ifdef CHECK_IDEAL_MWALK
5343      if(printout > 1)
5344      {
5345        idString(M,"//** Mwalk: M");
5346      }
5347//#endif
5348      rChangeCurrRing(baseRing);
5349      M1 =  idrMoveR(M, newRing,currRing);
5350      idDelete(&M);
5351      Gomega2 = idrMoveR(Gomega1, newRing,currRing);
5352      idDelete(&Gomega1);
5353      //PrintS("\n //** Mwalk: MLifttwoIdeal");
5354      F = MLifttwoIdeal(Gomega2, M1, G);
5355//#ifdef CHECK_IDEAL_MWALK
5356      if(printout > 2)
5357      {
5358        idString(F,"//** Mwalk: F");
5359      }
5360//#endif
5361      idDelete(&Gomega2);
5362      idDelete(&M1);
5363      rChangeCurrRing(newRing); // change the ring to newRing
5364      G = idrMoveR(F,baseRing,currRing);
5365      idDelete(&F);
5366      baseRing = currRing;
5367      si_opt_1 = save1; //set original options, e. g. option(RedSB)
5368      idSkipZeroes(G);
5369#ifdef TIME_TEST
5370      to = clock();
5371#endif
5372      //PrintS("\n //**Mwalk: Interreduce");
5373      //interreduce the Groebner basis <G> w.r.t. currRing
5374      //G = kInterRedCC(G,NULL);
5375#ifdef TIME_TEST
5376      tred = tred + clock() - to;
5377#endif
5378      idSkipZeroes(G);
5379      delete next_weight; */
5380      break;
5381    }
5382
5383    for(i=nV-1; i>=0; i--)
5384    {
5385      //(*tmp_weight)[i] = (*curr_weight)[i];
5386      (*curr_weight)[i] = (*next_weight)[i];
5387    }
5388    delete next_weight;
5389  }
5390  rChangeCurrRing(XXRing);
5391  ideal result = idrMoveR(G,baseRing,currRing);
5392  idDelete(&Go);
5393  idDelete(&G);
5394  //delete tmp_weight;
5395  delete ivNull;
5396  delete exivlp;
5397#ifndef BUCHBERGER_ALG
5398  delete last_omega;
5399#endif
5400#ifdef TIME_TEST
5401  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
5402  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
5403  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
5404#endif
5405  Print("\n//** Mwalk: Groebner Walk took %d steps.\n", nstep);
5406  return(result);
5407}
5408
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)
5412{
5413  BITSET save1 = si_opt_1; // save current options
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  }
5419  Set_Error(FALSE);
5420  Overflow_Error = FALSE;
5421  BOOLEAN endwalks = FALSE;
5422#ifdef TIME_TEST
5423  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
5424  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
5425  tinput = clock();
5426  clock_t tim;
5427#endif
5428  nstep=0;
5429  int i,polylength,nwalk;
5430  int nV = currRing->N;
5431
5432  ideal Gomega, M, F,FF, Gomega1, Gomega2, M1;
5433  ring newRing;
5434  ring targetRing;
5435  ring baseRing = currRing;
5436  ring XXRing = currRing;
5437  intvec* ivNull = new intvec(nV);
5438  intvec* curr_weight = new intvec(nV);
5439  intvec* target_weight = new intvec(nV);
5440  intvec* next_weight= new intvec(nV);
5441
5442  for(i=0; i<nV; i++)
5443  {
5444    (*curr_weight)[i] = (*orig_M)[i];
5445    (*target_weight)[i] = (*target_M)[i];
5446  }
5447#ifndef BUCHBERGER_ALG
5448  intvec* hilb_func;
5449   // to avoid (1,0,...,0) as the target vector
5450  intvec* last_omega = new intvec(nV);
5451  for(i=nV-1; i>0; i--)
5452  {
5453    (*last_omega)[i] = 1;
5454  }
5455  (*last_omega)[0] = 10000;
5456#endif
5457  rComplete(currRing);
5458#ifdef TIME_TEST
5459  to = clock();
5460#endif
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  }
5479  rChangeCurrRing(newRing);
5480  ideal G = MstdCC(idrMoveR(Go,baseRing,currRing));
5481  baseRing = currRing;
5482#ifdef TIME_TEST
5483  tostd = clock()-to;
5484#endif
5485
5486  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
5496  while(1)
5497  {
5498    nwalk ++;
5499    nstep ++;
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    }
5516#ifndef  BUCHBERGER_ALG
5517    if(isNolVector(curr_weight) == 0)
5518    {
5519      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
5520    }   
5521    else
5522    {
5523      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
5524    }
5525#endif
5526    if(nwalk == 1)
5527    {
5528      if(orig_M->length() == nV)
5529      {
5530        newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
5531      }
5532      else
5533      {
5534        newRing = VMatrDefault(orig_M);
5535      }
5536    }
5537    else
5538    {
5539     if(target_M->length() == nV)
5540     {
5541       newRing = VMrDefault(curr_weight); // define a new ring with ordering "(a(curr_weight),lp)
5542     }
5543     else
5544     {
5545       newRing = VMatrRefine(target_M,curr_weight);
5546     }
5547    }
5548    rChangeCurrRing(newRing);
5549    Gomega1 = idrMoveR(Gomega, baseRing,currRing);
5550    idDelete(&Gomega);
5551    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
5552#ifdef TIME_TEST
5553    to = clock();
5554#endif
5555#ifndef  BUCHBERGER_ALG
5556    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
5557    delete hilb_func;
5558#else
5559    M = kStd(Gomega1,NULL,testHomog,NULL,NULL,0,0,NULL);
5560#endif
5561#ifdef TIME_TEST
5562    tstd = tstd + clock() - to;
5563#endif
5564    idSkipZeroes(M);
5565//#ifdef CHECK_IDEAL_MWALK
5566    if(printout > 2)
5567    {
5568      idString(M, "//** Mrwalk: M");
5569    }
5570//#endif
5571    //change the ring to baseRing
5572    rChangeCurrRing(baseRing);
5573    M1 =  idrMoveR(M, newRing,currRing);
5574    idDelete(&M);
5575    Gomega2 = idrMoveR(Gomega1, newRing,currRing);
5576    idDelete(&Gomega1);
5577#ifdef TIME_TEST
5578    to = clock();
5579#endif
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
5582    F = MLifttwoIdeal(Gomega2, M1, G);
5583#ifdef TIME_TEST
5584    tlift = tlift + clock() - to;
5585#endif
5586//#ifdef CHECK_IDEAL_MWALK
5587    if(printout > 2)
5588    {
5589      idString(F, "//** Mrwalk: F");
5590    }
5591//#endif
5592    idDelete(&Gomega2);
5593    idDelete(&M1);
5594    rChangeCurrRing(newRing); // change the ring to newRing
5595    G = idrMoveR(F,baseRing,currRing);
5596    idDelete(&F);
5597    baseRing = currRing;
5598#ifdef TIME_TEST
5599    to = clock();
5600    tstd = tstd + clock() - to;
5601#endif
5602    idSkipZeroes(G);
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:
5625#ifdef TIME_TEST
5626    to = clock();
5627#endif
5628    next_weight = MwalkNextWeightCC(curr_weight,target_weight,G);
5629#ifdef TIME_TEST
5630    tnw = tnw + clock() - to;
5631#endif
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
5655      idDelete(&Gomega);
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    {
5668      baseRing = currRing;
5669      delete next_weight;
5670      break;
5671    }
5672
5673//#ifdef PRINT_VECTORS
5674    if(printout > 0)
5675    {
5676      MivString(curr_weight, target_weight, next_weight);
5677    }
5678//#endif
5679
5680    for(i=nV-1; i>=0; i--)
5681    {
5682      (*curr_weight)[i] = (*next_weight)[i];
5683    }
5684    delete next_weight;
5685  }
5686  baseRing = currRing;
5687  rChangeCurrRing(XXRing);
5688  ideal result = idrMoveR(G,baseRing,currRing);
5689  idDelete(&G);
5690  si_opt_1 = save1; //set original options, e. g. option(RedSB)
5691  delete ivNull;
5692#ifndef BUCHBERGER_ALG
5693  delete last_omega;
5694#endif
5695  Print("\n//** Mrwalk: Groebner Walk took %d steps.\n", nstep);
5696#ifdef TIME_TEST
5697  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
5698  //Print("\n// pSetm_Error = (%d)", ErrorCheck());
5699  //Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
5700#endif
5701  return(result);
5702}
5703
5704//unused
5705#if 0
5706ideal Mwalk_tst(ideal Go, intvec* curr_weight, intvec* target_weight)
5707{
5708  //clock_t tinput=clock();
5709  //idString(Go,"Ginp");
5710  int i, nV = currRing->N;
5711  int nwalk=0, endwalks=0;
5712
5713  ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G;
5714  // ideal G1; ring endRing;
5715  ring newRing, oldRing;
5716  intvec* ivNull = new intvec(nV);
5717  ring XXRing = currRing;
5718
5719  intvec* tmp_weight = new intvec(nV);
5720  for(i=nV-1; i>=0; i--)
5721  {
5722    (*tmp_weight)[i] = (*curr_weight)[i];
5723  }
5724  /* the monomial ordering of this current ring would be "dp" */
5725  G = MstdCC(Go);
5726#ifndef BUCHBERGER_ALG
5727  intvec* hilb_func;
5728#endif
5729  /* to avoid (1,0,...,0) as the target vector */
5730  intvec* last_omega = new intvec(nV);
5731  for(i=nV-1; i>0; i--)
5732    (*last_omega)[i] = 1;
5733  (*last_omega)[0] = 10000;
5734
5735  while(1)
5736  {
5737    nwalk ++;
5738    //Print("\n// Entering the %d-th step:", nwalk);
5739    //Print("\n// ring r[%d] = %s;", nwalk, rString(currRing));
5740    idString(G,"G");
5741    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
5742    Gomega = MwalkInitialForm(G, curr_weight);
5743    //ivString(curr_weight, "omega");
5744    idString(Gomega,"Gw");
5745
5746#ifndef  BUCHBERGER_ALG
5747    if(isNolVector(curr_weight) == 0)
5748      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
5749    else
5750      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
5751#endif // BUCHBERGER_ALG
5752
5753
5754    oldRing = currRing;
5755
5756    /* define a new ring that its ordering is "(a(curr_weight),lp) */
5757    VMrDefault(curr_weight);
5758    newRing = currRing;
5759
5760    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
5761
5762    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
5763#ifdef  BUCHBERGER_ALG
5764    M = MstdhomCC(Gomega1);
5765#else
5766    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
5767    delete hilb_func;
5768#endif // BUCHBERGER_ALG
5769
5770    idString(M,"M");
5771
5772      /* change the ring to oldRing */
5773    rChangeCurrRing(oldRing);
5774    M1 =  idrMoveR(M, newRing,currRing);
5775    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
5776
5777      /* compute a representation of the generators of submod (M)
5778         with respect to those of mod (Gomega).
5779         Gomega is a reduced Groebner basis w.r.t. the current ring */
5780    F = MLifttwoIdeal(Gomega2, M1, G);
5781    idDelete(&M1);
5782    idDelete(&Gomega2);
5783    idDelete(&G);
5784    idString(F,"F");
5785
5786    /* change the ring to newRing */
5787    rChangeCurrRing(newRing);
5788    F1 = idrMoveR(F, oldRing,currRing);
5789
5790    /* reduce the Groebner basis <G> w.r.t. new ring */
5791    G = kInterRedCC(F1, NULL);
5792    //idSkipZeroes(G);//done by kInterRed
5793    idDelete(&F1);
5794    idString(G,"G");
5795    if(endwalks == 1)
5796      break;
5797
5798    /* compute a next weight vector */
5799    intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G);
5800#ifdef PRINT_VECTORS
5801    MivString(curr_weight, target_weight, next_weight);
5802#endif
5803
5804    if(MivComp(next_weight, ivNull) == 1)
5805    {
5806      delete next_weight;
5807      break;
5808    }
5809    if(MivComp(next_weight, target_weight) == 1)
5810      endwalks = 1;
5811
5812    for(i=nV-1; i>=0; i--)
5813      (*tmp_weight)[i] = (*curr_weight)[i];
5814
5815    /* 06.11.01 to free the memory: NOT Changed!!*/
5816    for(i=nV-1; i>=0; i--)
5817      (*curr_weight)[i] = (*next_weight)[i];
5818    delete next_weight;
5819  }
5820  rChangeCurrRing(XXRing);
5821  G = idrMoveR(G, newRing,currRing);
5822
5823  delete tmp_weight;
5824  delete ivNull;
5825  PrintLn();
5826  return(G);
5827}
5828#endif
5829
5830/**************************************************************/
5831/*     Implementation of the perturbation walk algorithm      */
5832/**************************************************************/
5833/* If the perturbed target weight vector or an intermediate weight vector
5834   doesn't stay in the correct Groebner cone, we have only
5835   a reduced Groebner basis for the given ideal with respect to
5836   a monomial order which differs to the given order.
5837   Then we have to compute the wanted  reduced Groebner basis for it.
5838   For this, we can use
5839   1) the improved Buchberger algorithm or
5840   2) the changed perturbation walk algorithm with a decreased degree.
5841*/
5842// use kStd, if nP = 0, else call LastGB
5843ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,
5844             intvec* target_weight, int nP, int reduction, int printout)
5845{
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  }
5853  Set_Error(FALSE  );
5854  Overflow_Error = FALSE;
5855  //Print("// pSetm_Error = (%d)", ErrorCheck());
5856
5857  clock_t  tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
5858  xtextra=0;
5859  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
5860  tinput = clock();
5861
5862  clock_t tim;
5863
5864  nstep = 0;
5865  int i, ntwC=1, ntestw=1,  nV = currRing->N;
5866  BOOLEAN endwalks = FALSE;
5867
5868  ideal Gomega, M, F, FF, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
5869  ring newRing, oldRing, TargetRing;
5870  intvec* iv_M_dp;
5871  intvec* iv_M_lp;
5872  intvec* exivlp = Mivlp(nV);
5873  intvec* orig_target = target_weight;
5874  intvec* pert_target_vector = target_weight;
5875  intvec* ivNull = new intvec(nV);
5876  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
5877#ifndef BUCHBERGER_ALG
5878  intvec* hilb_func;
5879#endif
5880  intvec* next_weight;
5881
5882  // to avoid (1,0,...,0) as the target vector
5883  intvec* last_omega = new intvec(nV);
5884  for(i=nV-1; i>0; i--)
5885    (*last_omega)[i] = 1;
5886  (*last_omega)[0] = 10000;
5887
5888  ring XXRing = currRing;
5889
5890  to = clock();
5891  // perturbs the original vector
5892  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
5893  {
5894    G = MstdCC(Go);
5895    tostd = clock()-to;
5896    if(op_deg != 1){
5897      iv_M_dp = MivMatrixOrderdp(nV);
5898      //ivString(iv_M_dp, "iv_M_dp");
5899      curr_weight = MPertVectors(G, iv_M_dp, op_deg);
5900    }
5901  }
5902  else
5903  {
5904    //define ring order := (a(curr_weight),lp);
5905    if (rParameter(currRing) != NULL)
5906      DefRingPar(curr_weight);
5907    else
5908      rChangeCurrRing(VMrDefault(curr_weight));
5909
5910    G = idrMoveR(Go, XXRing,currRing);
5911    G = MstdCC(G);
5912    tostd = clock()-to;
5913    if(op_deg != 1){
5914      iv_M_dp = MivMatrixOrder(curr_weight);
5915      curr_weight = MPertVectors(G, iv_M_dp, op_deg);
5916    }
5917  }
5918  delete iv_dp;
5919  if(op_deg != 1) delete iv_M_dp;
5920
5921  ring HelpRing = currRing;
5922
5923  // perturbs the target weight vector
5924  if(tp_deg > 1 && tp_deg <= nV)
5925  {
5926    if (rParameter(currRing) != NULL)
5927      DefRingPar(target_weight);
5928    else
5929      rChangeCurrRing(VMrDefault(target_weight));
5930
5931    TargetRing = currRing;
5932    ssG = idrMoveR(G,HelpRing,currRing);
5933    if(MivSame(target_weight, exivlp) == 1)
5934    {
5935      iv_M_lp = MivMatrixOrderlp(nV);
5936      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
5937    }
5938    else
5939    {
5940      iv_M_lp = MivMatrixOrder(target_weight);
5941      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
5942    }
5943    delete iv_M_lp;
5944    pert_target_vector = target_weight;
5945    rChangeCurrRing(HelpRing);
5946    G = idrMoveR(ssG, TargetRing,currRing);
5947  }
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    }
5954  while(1)
5955  {
5956    nstep ++;
5957    to = clock();
5958    // compute an initial form ideal of <G> w.r.t. the weight vector
5959    // "curr_weight"
5960    Gomega = MwalkInitialForm(G, curr_weight);
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    }
6009
6010#ifdef ENDWALKS
6011    if(endwalks == TRUE)
6012    {
6013      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
6014      idElements(G, "G");
6015      headidString(G, "G");
6016    }
6017#endif
6018
6019    tif = tif + clock()-to;
6020
6021#ifndef  BUCHBERGER_ALG
6022    if(isNolVector(curr_weight) == 0)
6023      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
6024    else
6025      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
6026#endif // BUCHBERGER_ALG
6027
6028    oldRing = currRing;
6029
6030    // define a new ring with ordering "(a(curr_weight),lp)
6031    if (rParameter(currRing) != NULL)
6032      DefRingPar(curr_weight);
6033    else
6034      rChangeCurrRing(VMrDefault(curr_weight));
6035
6036    newRing = currRing;
6037    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
6038
6039#ifdef ENDWALKS
6040    if(endwalks==TRUE)
6041    {
6042      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
6043      idElements(Gomega1, "Gw");
6044      headidString(Gomega1, "headGw");
6045      PrintS("\n// compute a rGB of Gw:\n");
6046
6047#ifndef  BUCHBERGER_ALG
6048      ivString(hilb_func, "w");
6049#endif
6050    }
6051#endif
6052
6053    tim = clock();
6054    to = clock();
6055    // compute a reduced Groebner basis of <Gomega> w.r.t. "newRing"
6056#ifdef  BUCHBERGER_ALG
6057    M = MstdhomCC(Gomega1);
6058#else
6059    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
6060    delete hilb_func;
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){
6070      xtstd = xtstd+clock()-to;
6071#ifdef ENDWALKS
6072      Print("\n// time for the last std(Gw)  = %.2f sec\n",
6073            ((double) clock())/1000000 -((double)tim) /1000000);
6074#endif
6075    }
6076    else
6077      tstd=tstd+clock()-to;
6078
6079    // change the ring to oldRing
6080    rChangeCurrRing(oldRing);
6081    M1 =  idrMoveR(M, newRing,currRing);
6082    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
6083
6084    to=clock();
6085    /* compute a representation of the generators of submod (M)
6086       with respect to those of mod (Gomega).
6087       Gomega is a reduced Groebner basis w.r.t. the current ring */
6088    F = MLifttwoIdeal(Gomega2, M1, G);
6089    if(endwalks == FALSE)
6090      tlift = tlift+clock()-to;
6091    else
6092      xtlift=clock()-to;
6093
6094//#ifdef CHECK_IDEAL_MWALK
6095    if(printout > 2)
6096    {
6097      idString(F,"//** Mpwalk: F");
6098    }
6099//#endif
6100
6101    idDelete(&M1);
6102    idDelete(&Gomega2);
6103    idDelete(&G);
6104
6105    // change the ring to newRing
6106    rChangeCurrRing(newRing);
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:
6130    to=clock();
6131    // compute a next weight vector
6132    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
6133    tnw=tnw+clock()-to;
6134//#ifdef PRINT_VECTORS
6135    if(printout > 0)
6136    {
6137      MivString(curr_weight, target_weight, next_weight);
6138    }
6139//#endif
6140
6141    if(Overflow_Error == TRUE)
6142    {
6143      ntwC = 0;
6144      //ntestomega = 1;
6145      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
6146      //idElements(G, "G");
6147      delete next_weight;
6148      goto FINISH_160302;
6149    }
6150    if(MivComp(next_weight, ivNull) == 1){
6151      newRing = currRing;
6152      delete next_weight;
6153      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
6154      break;
6155    }
6156    if(MivComp(next_weight, target_weight) == 1)
6157      endwalks = TRUE;
6158
6159    for(i=nV-1; i>=0; i--)
6160      (*curr_weight)[i] = (*next_weight)[i];
6161
6162    delete next_weight;
6163  }//end of while-loop
6164
6165  if(tp_deg != 1)
6166  {
6167  FINISH_160302:
6168    if(MivSame(orig_target, exivlp) == 1)
6169      if (rParameter(currRing) != NULL)
6170        DefRingParlp();
6171      else
6172        VMrDefaultlp();
6173    else
6174      if (rParameter(currRing) != NULL)
6175        DefRingPar(orig_target);
6176      else
6177        rChangeCurrRing(VMrDefault(orig_target));
6178
6179    TargetRing=currRing;
6180    F1 = idrMoveR(G, newRing,currRing);
6181#ifdef CHECK_IDEAL
6182      headidString(G, "G");
6183#endif
6184
6185
6186    // check whether the pertubed target vector stays in the correct cone
6187    if(ntwC != 0){
6188      ntestw = test_w_in_ConeCC(F1, pert_target_vector);
6189    }
6190
6191    if( ntestw != 1 || ntwC == 0)
6192    {
6193      if(ntestw != 1 && printout >2)
6194      {
6195        ivString(pert_target_vector, "tau");
6196        PrintS("\n// ** perturbed target vector doesn't stay in cone!!");
6197        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
6198        //idElements(F1, "G");
6199      }
6200      // LastGB is "better" than the kStd subroutine
6201      to=clock();
6202      ideal eF1;
6203      if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1){
6204        // PrintS("\n// ** calls \"std\" to compute a GB");
6205        eF1 = MstdCC(F1);
6206        idDelete(&F1);
6207      }
6208      else {
6209        // PrintS("\n// ** calls \"LastGB\" to compute a GB");
6210        rChangeCurrRing(newRing);
6211        ideal F2 = idrMoveR(F1, TargetRing,currRing);
6212        eF1 = LastGB(F2, curr_weight, tp_deg-1);
6213        F2=NULL;
6214      }
6215      xtextra=clock()-to;
6216      ring exTargetRing = currRing;
6217
6218      rChangeCurrRing(XXRing);
6219      Eresult = idrMoveR(eF1, exTargetRing,currRing);
6220    }
6221    else{
6222      rChangeCurrRing(XXRing);
6223      Eresult = idrMoveR(F1, TargetRing,currRing);
6224    }
6225  }
6226  else {
6227    rChangeCurrRing(XXRing);
6228    Eresult = idrMoveR(G, newRing,currRing);
6229  }
6230  si_opt_1 = save1; //set original options, e. g. option(RedSB)
6231  delete ivNull;
6232  if(tp_deg != 1)
6233    delete target_weight;
6234
6235  if(op_deg != 1 )
6236    delete curr_weight;
6237
6238  delete exivlp;
6239  delete last_omega;
6240
6241#ifdef TIME_TEST
6242  TimeStringFractal(tinput, tostd, tif+xtif, tstd+xtstd,0, tlift+xtlift, tred+xtred,
6243             tnw+xtnw);
6244
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);
6777  return(Eresult);
6778}
6779
6780intvec* XivNull;
6781
6782// define a matrix (1 ... 1)
6783intvec* MMatrixone(int nV)
6784{
6785  int i,j;
6786  intvec* ivM = new intvec(nV*nV);
6787
6788  for(i=0; i<nV; i++)
6789    for(j=0; j<nV; j++)
6790    (*ivM)[i*nV + j] = 1;
6791
6792  return(ivM);
6793}
6794
6795int nnflow;
6796int Xcall;
6797int Xngleich;
6798
6799/***********************************************************************
6800 * Perturb the start weight vector at the top level, i.e. nlev = 1     *
6801 ***********************************************************************/
6802static ideal rec_fractal_call(ideal G, int nlev, intvec* ivtarget,
6803             int reduction, int printout)
6804{
6805  Overflow_Error =  FALSE;
6806  if(printout >0)
6807  {
6808    Print("\n\n// Entering the %d-th recursion:", nlev);
6809  }
6810  int i, nV = currRing->N;
6811  ring new_ring, testring;
6812  //ring extoRing;
6813  ideal Gomega, Gomega1, Gomega2, FF, F, F1, Gresult, Gresult1, G1, Gt;
6814  int nwalks = 0;
6815  intvec* Mwlp;
6816#ifndef BUCHBERGER_ALG
6817  intvec* hilb_func;
6818#endif
6819//  intvec* extXtau;
6820  intvec* next_vect;
6821  intvec* omega2 = new intvec(nV);
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  }
6829  //BOOLEAN isnewtarget = FALSE;
6830
6831  // to avoid (1,0,...,0) as the target vector (Hans)
6832  intvec* last_omega = new intvec(nV);
6833  for(i=nV-1; i>0; i--)
6834    (*last_omega)[i] = 1;
6835  (*last_omega)[0] = 10000;
6836
6837  intvec* omega = new intvec(nV);
6838  for(i