source: git/Singular/walk.cc @ d471a7

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