source: git/Singular/walk.cc @ d30a399

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