source: git/Singular/walk.cc @ 9614bb

spielwiese
Last change on this file since 9614bb was a601d5, checked in by Hans Schoenemann <hannes@…>, 11 years ago
chg: code simplification in walk
  • Property mode set to 100644
File size: 139.8 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 = 4;
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  /*weights: entries for 3 blocks: NULL Made:???*/
1854  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
1855  r->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
1856  for(i=0; i<nv; i++)
1857    r->wvhdl[0][i] = (*va)[i];
1858
1859  /* order: a,lp,C,0 */
1860  r->order = (int *) omAlloc(nb * sizeof(int *));
1861  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
1862  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
1863
1864  /* ringorder a for the first block: var 1..nv */
1865  r->order[0]  = ringorder_a;
1866  r->block0[0] = 1;
1867  r->block1[0] = nv;
1868
1869  /* ringorder lp for the second block: var 1..nv */
1870  r->order[1]  = ringorder_lp;
1871  r->block0[1] = 1;
1872  r->block1[1] = nv;
1873
1874  /* ringorder C for the third block */
1875  // it is very important within "idLift",
1876  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
1877  // therefore, nb  must be (nBlocks(currRing)  + 1)
1878  r->order[2]  = ringorder_C;
1879
1880  /* the last block: everything is 0 */
1881  r->order[3]  = 0;
1882
1883  /*polynomial ring*/
1884  r->OrdSgn    = 1;
1885
1886  /* complete ring intializations */
1887  rComplete(r);
1888
1889  rChangeCurrRing(r);
1890}
1891
1892/* 03.11.01 */
1893/* define and execute a new ring which order is  a lexicographic order */
1894static void VMrDefaultlp(void)
1895{
1896
1897  if ((currRing->ppNoether)!=NULL)
1898    pDelete(&(currRing->ppNoether));
1899
1900  if (((sLastPrinted.rtyp>BEGIN_RING) && (sLastPrinted.rtyp<END_RING)) ||
1901      ((sLastPrinted.rtyp==LIST_CMD)&&(lRingDependend((lists)sLastPrinted.data))))
1902
1903  {
1904    sLastPrinted.CleanUp();
1905  }
1906
1907  ring r = (ring) omAlloc0Bin(sip_sring_bin);
1908  int i, nv = currRing->N;
1909
1910  r->cf  = currRing->cf;
1911  r->N   = currRing->N;
1912  int nb = rBlocks(currRing) + 1;//31.10.01 (+1)
1913
1914  /*names*/
1915  char* Q; //30.10.01 to avoid the corrupted memory, NOT change!!
1916  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
1917  for(i=0; i<nv; i++)
1918  {
1919    Q = currRing->names[i];
1920    r->names[i]  = omStrDup(Q);
1921  }
1922
1923  /*weights: entries for 3 blocks: NULL Made:???*/
1924
1925  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
1926
1927  /* order: lp,C,0 */
1928  r->order = (int *) omAlloc(nb * sizeof(int *));
1929  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
1930  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
1931
1932  /* ringorder lp for the first block: var 1..nv */
1933  r->order[0]  = ringorder_lp;
1934  r->block0[0] = 1;
1935  r->block1[0] = nv;
1936
1937  /* ringorder C for the second block */
1938  r->order[1]  = ringorder_C;
1939
1940  /* the last block: everything is 0 */
1941  r->order[2]  = 0;
1942
1943  /*polynomial ring*/
1944  r->OrdSgn    = 1;
1945
1946  /* complete ring intializations */
1947  rComplete(r);
1948
1949  rChangeCurrRing(r);
1950}
1951
1952
1953/* define a ring with parameters und change to it */
1954/* DefRingPar and DefRingParlp corrupt still memory */
1955static void DefRingPar(intvec* va)
1956{
1957  int i, nv = currRing->N;
1958  int nb = rBlocks(currRing) + 1;
1959
1960  ring res=(ring)omAllocBin(sip_sring_bin);
1961
1962  memcpy(res,currRing,sizeof(ip_sring));
1963
1964  res->VarOffset = NULL;
1965  res->ref=0;
1966
1967  res->cf = currRing->cf; currRing->cf->ref++;
1968
1969//   if (currRing->cf->extRing!=NULL)
1970//     currRing->cf->extRing->ref++;
1971//
1972//   if (rParameter (currRing)!=NULL)
1973//   {
1974//     res->cf->extRing->qideal->m[0]=p_Copy(currRing->cf->extRing->qideal->m[0],currRing->cf->extRing);
1975//     int l=rPar(currRing);
1976//
1977//     res->cf->extRing->names=(char **)omAlloc(l*sizeof(char_ptr));
1978//
1979//     for(i=l-1;i>=0;i--)
1980//       rParameter (res)[i]=omStrDup(rParameter (currRing)[i]);
1981//   }
1982
1983  /*weights: entries for 3 blocks: NULL Made:???*/
1984  res->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
1985  res->wvhdl[0] = (int*) omAlloc(nv*sizeof(int));
1986  for(i=0; i<nv; i++)
1987    res->wvhdl[0][i] = (*va)[i];
1988
1989  /* order: a,lp,C,0 */
1990  res->order = (int *) omAlloc(nb * sizeof(int *));
1991  res->block0 = (int *)omAlloc0(nb * sizeof(int *));
1992  res->block1 = (int *)omAlloc0(nb * sizeof(int *));
1993
1994  /* ringorder a for the first block: var 1..nv */
1995  res->order[0]  = ringorder_a;
1996  res->block0[0] = 1;
1997  res->block1[0] = nv;
1998
1999  /* ringorder lp for the second block: var 1..nv */
2000  res->order[1]  = ringorder_lp;
2001  res->block0[1] = 1;
2002  res->block1[1] = nv;
2003
2004  /* ringorder C for the third block */
2005  // it is very important within "idLift",
2006  // especially, by ring syz_ring=rCurrRingAssure_SyzComp();
2007  // therefore, nb  must be (nBlocks(currRing)  + 1)
2008  res->order[2]  = ringorder_C;
2009
2010  /* the last block: everything is 0 */
2011  res->order[3]  = 0;
2012
2013  /*polynomial ring*/
2014  res->OrdSgn    = 1;
2015
2016
2017  res->names   = (char **)omAlloc0(nv * sizeof(char_ptr));
2018  for (i=nv-1; i>=0; i--)
2019    res->names[i] = omStrDup(currRing->names[i]);
2020
2021  /* complete ring intializations */
2022   rComplete(res);
2023
2024
2025  // clean up history
2026  if (sLastPrinted.RingDependend())
2027  {
2028    sLastPrinted.CleanUp();
2029  }
2030
2031
2032  /* execute the created ring */
2033  rChangeCurrRing(res);
2034}
2035
2036
2037static void DefRingParlp(void)
2038{
2039  int i, nv = currRing->N;
2040
2041  ring r=(ring)omAllocBin(sip_sring_bin);
2042
2043  memcpy(r,currRing,sizeof(ip_sring));
2044
2045  r->VarOffset = NULL;
2046  r->ref=0;
2047
2048  r->cf = currRing->cf; currRing->cf->ref++;
2049
2050//   if (currRing->cf->extRing!=NULL)
2051//     currRing->cf->extRing->ref++;
2052//
2053//   if (rParameter (currRing)!=NULL)
2054//   {
2055//     r->cf->extRing->qideal->m[0]=p_Copy(currRing->cf->extRing->qideal->m[0], currRing->cf->extRing);
2056//     int l=rPar(currRing);
2057//     r->cf->extRing->names=(char **)omAlloc(l*sizeof(char_ptr));
2058//
2059//     for(i=l-1;i>=0;i--)
2060//       rParameter(r)[i]=omStrDup(rParameter (currRing)[i]);
2061//   }
2062
2063
2064  r->cf  = currRing->cf;
2065  r->N   = currRing->N;
2066  int nb = rBlocks(currRing) + 1;//31.10.01 (+1)
2067
2068  /*names*/
2069  char* Q;
2070  r->names = (char **) omAlloc0(nv * sizeof(char_ptr));
2071  for(i=nv-1; i>=0; i--)
2072  {
2073    Q = currRing->names[i];
2074    r->names[i]  = omStrDup(Q);
2075  }
2076
2077  /*weights: entries for 3 blocks: NULL Made:???*/
2078
2079  r->wvhdl = (int **)omAlloc0(nb * sizeof(int_ptr));
2080
2081  /* order: lp,C,0 */
2082  r->order = (int *) omAlloc(nb * sizeof(int *));
2083  r->block0 = (int *)omAlloc0(nb * sizeof(int *));
2084  r->block1 = (int *)omAlloc0(nb * sizeof(int *));
2085
2086  /* ringorder lp for the first block: var 1..nv */
2087  r->order[0]  = ringorder_lp;
2088  r->block0[0] = 1;
2089  r->block1[0] = nv;
2090
2091  /* ringorder C for the second block */
2092  r->order[1]  = ringorder_C;
2093
2094  /* the last block: everything is 0 */
2095  r->order[2]  = 0;
2096
2097  /*polynomial ring*/
2098  r->OrdSgn    = 1;
2099
2100
2101//   if (rParameter(currRing)!=NULL)
2102//   {
2103//     r->cf->extRing->qideal->m[0]=p_Copy(currRing->cf->extRing->qideal->m[0], currRing->cf->extRing);
2104//     int l=rPar(currRing);
2105//     r->cf->extRing->names=(char **)omAlloc(l*sizeof(char_ptr));
2106//
2107//     for(i=l-1;i>=0;i--)
2108//       rParameter(r)[i]=omStrDup(rParameter(currRing)[i]);
2109//   }
2110
2111  /* complete ring intializations */
2112  rComplete(r);
2113
2114  // clean up history
2115  if (sLastPrinted.RingDependend())
2116  {
2117    sLastPrinted.CleanUp();
2118  }
2119
2120  /* execute the created ring */
2121  rChangeCurrRing(r);
2122}
2123
2124#if 0 /*unused*/
2125/* check wheather one or more components of a vector are zero */
2126static int isNolVector(intvec* hilb)
2127{
2128  int i;
2129  for(i=hilb->length()-1; i>=0; i--)
2130    if((* hilb)[i]==0)
2131      return 1;
2132
2133  return 0;
2134}
2135#endif
2136
2137
2138/******************************  Februar 2002  ****************************
2139 * G is a Groebner basis w.r.t. (a(curr_weight),lp) and                   *
2140 * we compute a GB of <G> w.r.t. the lex. order by the perturbation walk  *
2141 * its perturbation degree is tp_deg                                      *
2142 * We call the following subfunction LastGB, if                           *
2143 *     the computed intermediate weight vector or                         *
2144 *     the perturbed target weight vector                                 *
2145 * does NOT in the correct cone.                                          *
2146 **************************************************************************/
2147
2148static ideal LastGB(ideal G, intvec* curr_weight,int tp_deg)
2149{
2150  BOOLEAN nError = Overflow_Error;
2151  Overflow_Error = FALSE;
2152
2153  int i, nV = currRing->N;
2154  int nwalk=0, endwalks=0, nnwinC=1;
2155  int nlast = 0;
2156  ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG;
2157  ring newRing, oldRing, TargetRing;
2158  intvec* iv_M_lp;
2159  intvec* target_weight;
2160  intvec* iv_lp = Mivlp(nV); //define (1,0,...,0)
2161  intvec* pert_target_vector;
2162  intvec* ivNull = new intvec(nV);
2163  intvec* extra_curr_weight = new intvec(nV);
2164  intvec* hilb_func;
2165  intvec* next_weight;
2166
2167  /* to avoid (1,0,...,0) as the target vector */
2168  intvec* last_omega = new intvec(nV);
2169  for(i=nV-1; i>0; i--)
2170    (*last_omega)[i] = 1;
2171  (*last_omega)[0] = 10000;
2172
2173  ring EXXRing = currRing;
2174
2175  /* compute a pertubed weight vector of the target weight vector */
2176  if(tp_deg > 1 && tp_deg <= nV)
2177  {
2178    //..25.03.03 VMrDefaultlp();//    VMrDefault(target_weight);
2179    if (rParameter (currRing) != NULL)
2180      DefRingParlp();
2181    else
2182      VMrDefaultlp();
2183
2184    TargetRing = currRing;
2185    ssG = idrMoveR(G,EXXRing,currRing);
2186    iv_M_lp = MivMatrixOrderlp(nV);
2187    //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
2188    target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
2189    delete iv_M_lp;
2190    pert_target_vector = target_weight;
2191
2192    rChangeCurrRing(EXXRing);
2193    G = idrMoveR(ssG, TargetRing,currRing);
2194  }
2195  else
2196    target_weight = Mivlp(nV);
2197
2198  //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
2199
2200  while(1)
2201  {
2202    nwalk++;
2203    nstep++;
2204    to=clock();
2205    /* compute a next weight vector */
2206    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
2207    xtnw=xtnw+clock()-to;
2208#ifdef PRINT_VECTORS
2209    MivString(curr_weight, target_weight, next_weight);
2210#endif
2211
2212    if(Overflow_Error == TRUE){
2213      newRing = currRing;
2214      nnwinC = 0;
2215      if(tp_deg == 1)
2216        nlast = 1;
2217      delete next_weight;
2218
2219      //idElements(G, "G");
2220      //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
2221
2222      break;
2223    }
2224
2225    if(MivComp(next_weight, ivNull) == 1){
2226      //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
2227      newRing = currRing;
2228      delete next_weight;
2229      break;
2230    }
2231
2232    if(MivComp(next_weight, target_weight) == 1)
2233      endwalks = 1;
2234
2235    for(i=nV-1; i>=0; i--)
2236        (*extra_curr_weight)[i] = (*curr_weight)[i];
2237
2238    /* 06.11.01 NOT Changed */
2239    for(i=nV-1; i>=0; i--)
2240      (*curr_weight)[i] = (*next_weight)[i];
2241
2242    oldRing = currRing;
2243    to=clock();
2244    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
2245    Gomega = MwalkInitialForm(G, curr_weight);
2246    xtif=xtif+clock()-to;
2247
2248#ifdef ENDWALKS
2249    if(endwalks == 1){
2250      Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
2251      idElements(Gomega, "Gw");
2252      headidString(Gomega, "Gw");
2253    }
2254#endif
2255
2256#ifndef  BUCHBERGER_ALG
2257    if(isNolVector(curr_weight) == 0)
2258      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
2259    else
2260      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
2261#endif // BUCHBERGER_ALG
2262
2263    /* define a new ring that its ordering is "(a(curr_weight),lp) */
2264    //..25.03.03 VMrDefault(curr_weight);
2265    if (rParameter (currRing) != NULL)
2266      DefRingPar(curr_weight);
2267    else
2268      VMrDefault(curr_weight);
2269
2270    newRing = currRing;
2271    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
2272
2273    to=clock();
2274    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
2275#ifdef  BUCHBERGER_ALG
2276    M = MstdhomCC(Gomega1);
2277#else
2278    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
2279    delete hilb_func;
2280#endif // BUCHBERGER_ALG
2281    xtstd=xtstd+clock()-to;
2282    /* change the ring to oldRing */
2283    rChangeCurrRing(oldRing);
2284    M1 =  idrMoveR(M, newRing,currRing);
2285    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
2286
2287    to=clock();
2288    /* compute a reduced Groebner basis of <G> w.r.t. "newRing" */
2289    F = MLifttwoIdeal(Gomega2, M1, G);
2290    xtlift=xtlift+clock()-to;
2291
2292    idDelete(&M1);
2293    idDelete(&G);
2294
2295    /* change the ring to newRing */
2296    rChangeCurrRing(newRing);
2297    F1 = idrMoveR(F, oldRing,currRing);
2298
2299    to=clock();
2300    /* reduce the Groebner basis <G> w.r.t. new ring */
2301    G = kInterRedCC(F1, NULL);
2302    xtred=xtred+clock()-to;
2303    idDelete(&F1);
2304
2305    if(endwalks == 1){
2306      //Print("\n// ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
2307      break;
2308    }
2309
2310    delete next_weight;
2311  }//while
2312
2313  delete ivNull;
2314
2315  if(tp_deg != 1)
2316  {
2317    //..25.03.03 VMrDefaultlp();//define and execute the ring "lp"
2318    if (rParameter (currRing) != NULL)
2319      DefRingParlp();
2320    else
2321      VMrDefaultlp();
2322
2323    F1 = idrMoveR(G, newRing,currRing);
2324
2325    if(nnwinC == 0 || test_w_in_ConeCC(F1, pert_target_vector) != 1)
2326    {
2327      oldRing = currRing;
2328      rChangeCurrRing(newRing);
2329      G = idrMoveR(F1, oldRing,currRing);
2330      Print("\n// takes %d steps and calls the recursion of level %d:",
2331             nwalk, tp_deg-1);
2332
2333      F1 = LastGB(G,curr_weight, tp_deg-1);
2334    }
2335
2336    TargetRing = currRing;
2337    rChangeCurrRing(EXXRing);
2338    result = idrMoveR(F1, TargetRing,currRing);
2339  }
2340  else
2341  {
2342    if(nlast == 1)
2343    {
2344      //OMEGA_OVERFLOW_LASTGB:
2345      /*
2346      if(MivSame(curr_weight, iv_lp) == 1)
2347        if (rParameter(currRing) != NULL)
2348          DefRingParlp();
2349        else
2350          VMrDefaultlp();
2351      else
2352        if (rParameter(currRing) != NULL)
2353          DefRingPar(curr_weight);
2354        else
2355          VMrDefault(curr_weight);
2356      */
2357
2358        //..25.03.03 VMrDefaultlp();//define and execute the ring "lp"
2359        if (rParameter (currRing) != NULL)
2360          DefRingParlp();
2361        else
2362          VMrDefaultlp();
2363
2364
2365      F1 = idrMoveR(G, newRing,currRing);
2366      //Print("\n// Apply \"std\" in ring r%d_%d = %s;\n", tp_deg, nwalk, rString(currRing));
2367
2368      G = MstdCC(F1);
2369      idDelete(&F1);
2370      newRing = currRing;
2371    }
2372
2373    rChangeCurrRing(EXXRing);
2374    result = idrMoveR(G, newRing,currRing);
2375  }
2376  delete target_weight;
2377  delete last_omega;
2378  delete iv_lp;
2379
2380  if(Overflow_Error == FALSE)
2381    Overflow_Error = nError;
2382
2383  return(result);
2384}
2385
2386
2387/* check whether a polynomial of G has least 3 monomials */
2388static int lengthpoly(ideal G)
2389{
2390  int i;
2391  for(i=IDELEMS(G)-1; i>=0; i--)
2392#if 0
2393    if(pLength(G->m[i])>2)
2394      return 1;
2395#else
2396    if((G->m[i]!=NULL) /* len >=0 */
2397       && (G->m[i]->next!=NULL) /* len >=1 */
2398       && (G->m[i]->next->next!=NULL) /* len >=2 */
2399       && (G->m[i]->next->next->next!=NULL) /* len >=3 */
2400      //&& (G->m[i]->next->next->next->next!=NULL) /* len >=4 */
2401       ) return 1;
2402#endif
2403  return 0;
2404}
2405
2406/* check whether a polynomial of G has least 2 monomials */
2407static int islengthpoly2(ideal G)
2408{
2409  int i;
2410  for(i=IDELEMS(G)-1; i>=0; i--)
2411    if((G->m[i]!=NULL) /* len >=0 */
2412       && (G->m[i]->next!=NULL) /* len >=1 */
2413       && (G->m[i]->next->next!=NULL)) /* len >=2 */
2414      return 1;
2415
2416  return 0;
2417}
2418
2419
2420
2421/* Implementation of the improved Groebner walk algorithm which is written
2422   by Quoc-Nam Tran (2000).
2423   One perturbs the original target weight vector, only if
2424   the next intermediate weight vector is equal to the current target weight
2425   vector. This must be repeated until the wanted reduced Groebner basis
2426   to reach.
2427   If the numbers of variables is big enough, the representation of the origin
2428   weight vector may be very big. Therefore, it is possible the intermediate
2429   weight vector doesn't stay in the correct Groebner cone.
2430   In this case we have just a reduced Groebner basis of the given ideal
2431   with respect to another monomial order. Then we have to compute
2432   a wanted reduced Groebner basis of it with respect to the given order.
2433   At the following subroutine we use the improved Buchberger algorithm or
2434   the changed perturbation walk algorithm with a decrased degree.
2435 */
2436
2437/*2
2438* return the initial term of an ideal
2439*/
2440static ideal idHeadCC(ideal h)
2441{
2442  int i, nH =IDELEMS(h);
2443
2444  ideal m = idInit(nH,h->rank);
2445
2446  for (i=nH-1;i>=0; i--)
2447  {
2448    if (h->m[i]!=NULL)
2449      m->m[i]=pHead(h->m[i]);
2450  }
2451  return m;
2452}
2453
2454/* check whether two head-ideals are the same */
2455static inline int test_G_GB_walk(ideal H0, ideal H1)
2456{
2457  int i, nG = IDELEMS(H0);
2458
2459  if(nG != IDELEMS(H1))
2460    return 0;
2461
2462  for(i=nG-1; i>=0; i--)
2463  {
2464#if 0
2465    poly t;
2466    if((t=pSub(pCopy(H0->m[i]), pCopy(H1->m[i]))) != NULL)
2467    {
2468      pDelete(&t);
2469      return 0;
2470    }
2471    pDelete(&t);
2472#else
2473    if(!pEqualPolys(H0->m[i],H1->m[i]))
2474      return 0;
2475#endif
2476  }
2477
2478  return 1;
2479}
2480
2481#if 0 /*unused*/
2482/* 19.11.01 */
2483/* find the maximal total degree of polynomials in G */
2484static int Trandegreebound(ideal G)
2485{
2486  int i, nG = IDELEMS(G);
2487  int np=1, nV = currRing->N;
2488  int degtmp, result = 0;
2489  intvec* ivUnit = Mivdp(nV);
2490
2491  for(i=nG-1; i>=0; i--)
2492  {
2493    /* find the maximal total degree of the polynomial G[i] */
2494    degtmp = MwalkWeightDegree(G->m[i], ivUnit);
2495    if(degtmp > result)
2496      result = degtmp;
2497  }
2498  delete ivUnit;
2499  return result;
2500}
2501#endif
2502
2503/* perturb the weight vector iva w.r.t. the ideal G.
2504   the monomial order of the current ring is the w_1 weight lex. order.
2505   define w := d^(n-1)w_1+ d^(n-2)w_2, ...+ dw_(n-1)+ w_n
2506   where d := 1 + max{totdeg(g):g in G}*m, or
2507   d := (2*maxdeg*maxdeg + (nV+1)*maxdeg)*m;
2508*/
2509
2510#if 0 /*unused*/
2511//GMP
2512static intvec* TranPertVector(ideal G, intvec* iva)
2513{
2514  BOOLEAN nError = Overflow_Error;
2515  Overflow_Error = FALSE;
2516
2517  int i, j,  nG = IDELEMS(G);
2518  int nV = currRing->N;
2519
2520  // define the sequence which expresses the current monomial ordering
2521  // w_1 = iva; w_2 = (1,0,..,0); w_n = (0,...,0,1,0)
2522  intvec* ivMat = MivMatrixOrder(iva);
2523
2524  int  mtmp, m=(*iva)[0];
2525
2526  for(i=ivMat->length(); i>=0; i--)
2527  {
2528    mtmp = (*ivMat)[i];
2529
2530    if(mtmp <0) mtmp = -mtmp;
2531
2532    if(mtmp > m)
2533      m = mtmp;
2534  }
2535
2536  /* define the maximal total degree of polynomials of G */
2537  mpz_t ndeg;
2538  mpz_init(ndeg);
2539
2540 // 12 Juli 03
2541#ifndef UPPER_BOUND
2542  mpz_set_si(ndeg, Trandegreebound(G)+1);
2543#else
2544  mpz_t ztmp;
2545  mpz_init(ztmp);
2546
2547  mpz_t maxdeg;
2548  mpz_init_set_si(maxdeg, Trandegreebound(G));
2549
2550  //ndeg = (2*maxdeg*maxdeg + (nV+1)*maxdeg)*m;//Kalkbrenner (1999)
2551  mpz_pow_ui(ztmp, maxdeg, 2);
2552  mpz_mul_ui(ztmp, ztmp, 2);
2553  mpz_mul_ui(maxdeg, maxdeg, nV+1);
2554  mpz_add(ndeg, ztmp, maxdeg);
2555  mpz_mul_ui(ndeg, ndeg, m);
2556
2557  //PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m ");
2558  //Print("\n//         where d = %d, n = %d and bound = %d", maxdeg, nV, ndeg);
2559#endif //UPPER_BOUND
2560
2561  /* 29.08.03*/
2562#ifdef INVEPS_SMALL_IN_TRAN
2563 if(mpz_cmp_ui(ndeg, nV)>0 && nV > 3)
2564    mpz_cdiv_q_ui(ndeg, ndeg, nV);
2565
2566 //PrintS("\n// choose the \"small\" inverse epsilon:");
2567 //mpz_out_str(stdout, 10, ndeg);
2568#endif
2569  mpz_t deg_tmp;
2570  mpz_init_set(deg_tmp, ndeg);
2571
2572  mpz_t *ivres=( mpz_t *) omAlloc(nV*sizeof(mpz_t));
2573  mpz_init_set_si(ivres[nV-1],1);
2574
2575  for(i=nV-2; i>=0; i--)
2576  {
2577    mpz_init_set(ivres[i], deg_tmp);
2578    mpz_mul(deg_tmp, deg_tmp, ndeg);
2579  }
2580
2581  mpz_t *ivtmp=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
2582  for(i=0; i<nV; i++)
2583    mpz_init(ivtmp[i]);
2584
2585  mpz_t sing_int;
2586  mpz_init_set_ui(sing_int,  2147483647);
2587
2588  intvec* repr_vector = new intvec(nV);
2589
2590  /* define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n */
2591  for(i=0; i<nV; i++)
2592    for(j=0; j<nV; j++)
2593    {
2594      if( (*ivMat)[i*nV+j] >= 0 )
2595        mpz_mul_ui(ivres[i], ivres[i], (*ivMat)[i*nV+j]);
2596      else
2597      {
2598        mpz_mul_ui(ivres[i], ivres[i], -(*ivMat)[i*nV+j]);
2599        mpz_neg(ivres[i], ivres[i]);
2600      }
2601      mpz_add(ivtmp[j], ivtmp[j], ivres[i]);
2602    }
2603
2604  delete ivMat;
2605
2606  int ntrue=0;
2607  for(i=0; i<nV; i++)
2608  {
2609    (*repr_vector)[i] = mpz_get_si(ivtmp[i]);
2610    if(mpz_cmp(ivtmp[i], sing_int)>=0)
2611    {
2612      ntrue++;
2613      if(Overflow_Error == FALSE)
2614      {
2615        Overflow_Error = TRUE;
2616
2617        PrintS("\n// ** OVERFLOW in \"Repr.Vector\": ");
2618        mpz_out_str( stdout, 10, ivtmp[i]);
2619        PrintS(" is greater than 2147483647 (max. integer representation)");
2620        Print("\n//  So vector[%d] := %d is wrong!!\n",i+1,(*repr_vector)[i]);
2621      }
2622    }
2623  }
2624  if(Overflow_Error == TRUE)
2625  {
2626    ivString(repr_vector, "repvector");
2627    Print("\n// %d element(s) of it are overflow!!", ntrue);
2628  }
2629
2630  if(Overflow_Error == FALSE)
2631    Overflow_Error=nError;
2632
2633  omFree(ivres);
2634  omFree(ivtmp);
2635  return repr_vector;
2636}
2637#endif
2638
2639
2640
2641#if 0 /*unused*/
2642static intvec* TranPertVector_lp(ideal G)
2643{
2644  BOOLEAN nError = Overflow_Error;
2645  Overflow_Error = FALSE;
2646
2647  int i, j,  nG = IDELEMS(G);
2648  int nV = currRing->N;
2649
2650  /* define the maximal total degree of polynomials of G */
2651  mpz_t ndeg;
2652  mpz_init(ndeg);
2653
2654 // 12 Juli 03
2655#ifndef UPPER_BOUND
2656  mpz_set_si(ndeg, Trandegreebound(G)+1);
2657#else
2658  mpz_t ztmp;
2659  mpz_init(ztmp);
2660
2661  mpz_t maxdeg;
2662  mpz_init_set_si(maxdeg, Trandegreebound(G));
2663
2664  //ndeg = (2*maxdeg*maxdeg + (nV+1)*maxdeg);//Kalkbrenner (1999)
2665  mpz_pow_ui(ztmp, maxdeg, 2);
2666  mpz_mul_ui(ztmp, ztmp, 2);
2667  mpz_mul_ui(maxdeg, maxdeg, nV+1);
2668  mpz_add(ndeg, ztmp, maxdeg);
2669  /*
2670    PrintS("\n// with the new upper degree bound (2d^2+(n+1)d)*m ");
2671    Print("\n//         where d = %d, n = %d and bound = %d",
2672    mpz_get_si(maxdeg), nV, mpz_get_si(ndeg));
2673  */
2674#endif //UPPER_BOUND
2675
2676#ifdef INVEPS_SMALL_IN_TRAN
2677 if(mpz_cmp_ui(ndeg, nV)>0 && nV > 3)
2678    mpz_cdiv_q_ui(ndeg, ndeg, nV);
2679
2680 //PrintS("\n// choose the \"small\" inverse epsilon:");
2681 // mpz_out_str(stdout, 10, ndeg);
2682#endif
2683
2684  mpz_t deg_tmp;
2685  mpz_init_set(deg_tmp, ndeg);
2686
2687  mpz_t *ivres=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
2688  mpz_init_set_si(ivres[nV-1], 1);
2689
2690  for(i=nV-2; i>=0; i--)
2691  {
2692    mpz_init_set(ivres[i], deg_tmp);
2693    mpz_mul(deg_tmp, deg_tmp, ndeg);
2694  }
2695
2696  mpz_t sing_int;
2697  mpz_init_set_ui(sing_int,  2147483647);
2698
2699  intvec* repr_vector = new intvec(nV);
2700  int ntrue=0;
2701  for(i=0; i<nV; i++)
2702  {
2703    (*repr_vector)[i] = mpz_get_si(ivres[i]);
2704
2705    if(mpz_cmp(ivres[i], sing_int)>=0)
2706    {
2707      ntrue++;
2708      if(Overflow_Error == FALSE)
2709      {
2710        Overflow_Error = TRUE;
2711        PrintS("\n// ** OVERFLOW in \"Repr.Vector\": ");
2712        mpz_out_str( stdout, 10, ivres[i]);
2713        PrintS(" is greater than 2147483647 (max. integer representation)");
2714        Print("\n//  So vector[%d] := %d is wrong!!\n",i+1,(*repr_vector)[i]);
2715      }
2716    }
2717  }
2718  if(Overflow_Error == TRUE)
2719  {
2720    ivString(repr_vector, "repvector");
2721    Print("\n// %d element(s) of it are overflow!!", ntrue);
2722  }
2723  if(Overflow_Error == FALSE)
2724    Overflow_Error = nError;
2725
2726  omFree(ivres);
2727  return repr_vector;
2728}
2729#endif
2730
2731
2732#if 0 /*unused*/
2733//GMP
2734static intvec* RepresentationMatrix_Dp(ideal G, intvec* M)
2735{
2736  BOOLEAN nError = Overflow_Error;
2737  Overflow_Error = FALSE;
2738
2739  int i, j;
2740  int nV = currRing->N;
2741
2742  intvec* ivUnit = Mivdp(nV);
2743  int degtmp, maxdeg = 0;
2744
2745  for(i=IDELEMS(G)-1; i>=0; i--)
2746  {
2747    /* find the maximal total degree of the polynomial G[i] */
2748    degtmp = MwalkWeightDegree(G->m[i], ivUnit);
2749    if(degtmp > maxdeg)
2750      maxdeg = degtmp;
2751  }
2752
2753  mpz_t ztmp;
2754  mpz_init_set_si(ztmp, maxdeg);
2755  mpz_t *ivres=(mpz_t *)omAlloc(nV*sizeof(mpz_t));
2756  mpz_init_set_si(ivres[nV-1], 1); // (*ivres)[nV-1] = 1;
2757
2758  for(i=nV-2; i>=0; i--)
2759  {
2760    mpz_init_set(ivres[i], ztmp); //(*ivres)[i] = ztmp;
2761    mpz_mul_ui(ztmp, ztmp, maxdeg); //ztmp *=maxdeg;
2762  }
2763
2764  mpz_t *ivtmp=(mpz_t*)omAlloc(nV*sizeof(mpz_t));
2765  for(i=0; i<nV; i++)
2766    mpz_init(ivtmp[i]);
2767
2768  /* define ivtmp := ndeg^(n-1).w_1 + ndeg^(n-2).w_2 + ... + w_n */
2769  for(i=0; i<nV; i++)
2770    for(j=0; j<nV; j++)
2771    {
2772      if((*M)[i*nV+j] < 0)
2773      {
2774        mpz_mul_ui(ztmp, ivres[i], -(*M)[i*nV+j]);
2775        mpz_neg(ztmp, ztmp);
2776      }
2777      else
2778        mpz_mul_ui(ztmp, ivres[i], (*M)[i*nV+j]);
2779
2780      mpz_add(ivtmp[j], ivtmp[j], ztmp);
2781    }
2782
2783  mpz_t sing_int;
2784  mpz_init_set_ui(sing_int,  2147483647);
2785
2786  int ntrue=0;
2787  intvec* repvector = new intvec(nV);
2788  for(i=0; i<nV; i++)
2789  {
2790    (*repvector)[i] = mpz_get_si(ivtmp[i]);
2791    if(mpz_cmp(ivtmp[i], sing_int)>0)
2792    {
2793      ntrue++;
2794      if(Overflow_Error == FALSE)
2795      {
2796        Overflow_Error = TRUE;
2797        PrintS("\n// ** OVERFLOW in \"Repr.Matrix\": ");
2798        mpz_out_str( stdout, 10, ivtmp[i]);
2799        PrintS(" is greater than 2147483647 (max. integer representation)");
2800        Print("\n//  So vector[%d] := %d is wrong!!\n",i+1,(*repvector)[i]);
2801      }
2802    }
2803  }
2804  if(Overflow_Error == TRUE)
2805  {
2806    ivString(repvector, "repvector");
2807    Print("\n// %d element(s) of it are overflow!!", ntrue);
2808  }
2809
2810  if(Overflow_Error == FALSE)
2811    Overflow_Error = nError;
2812
2813  mpz_clear(ztmp);
2814  omFree(ivtmp);
2815  omFree(ivres);
2816  return repvector;
2817}
2818#endif
2819
2820
2821
2822
2823
2824/* The following subroutine is the implementation of our first improved
2825   Groebner walk algorithm, i.e. the first altervative algorithm.
2826   First we use the Grobner walk algorithm and then we call the changed
2827   perturbation walk algorithm with decreased degree, if an intermediate
2828   weight vector is equal to the current target weight vector.
2829   This call will be only repeated until we get the wanted reduced Groebner
2830   basis or n times, where n is the numbers of variables.
2831*/
2832
2833// 19 Juni 2003
2834#if 0 /* unused*/
2835static int testnegintvec(intvec* v)
2836{
2837  int n = v->length();
2838  int i;
2839  for(i=0; i<n; i++){
2840    if((*v)[i]<0)
2841      return(1);
2842  }
2843  return(0);
2844}
2845#endif
2846
2847
2848/* 7 Februar 2002 */
2849/* npwinc = 0, if curr_weight doesn't stay in the correct Groebner cone */
2850static ideal Rec_LastGB(ideal G, intvec* curr_weight,
2851                        intvec* orig_target_weight, int tp_deg, int npwinc)
2852{
2853  BOOLEAN nError = Overflow_Error;
2854  Overflow_Error = FALSE;
2855  BOOLEAN nOverflow_Error = FALSE;
2856
2857  clock_t tproc=0;
2858  clock_t tinput = clock();
2859
2860  int i,  nV = currRing->N;
2861  int nwalk=0, endwalks=0, nnwinC=1;
2862  int nlast = 0;
2863  ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG;
2864  ring newRing, oldRing, TargetRing;
2865  intvec* iv_M_lp;
2866  intvec* target_weight;
2867  intvec* ivNull = new intvec(nV); //define (0,...,0)
2868  ring EXXRing = currRing;
2869  int NEG=0; //19 juni 03
2870  intvec* extra_curr_weight = new intvec(nV);
2871  intvec* next_weight;
2872
2873  //08 Juli 03
2874  intvec* hilb_func;
2875  /* to avoid (1,0,...,0) as the target vector */
2876  intvec* last_omega = new intvec(nV);
2877  for(i=nV-1; i>0; i--)
2878    (*last_omega)[i] = 1;
2879  (*last_omega)[0] = 10000;
2880
2881  BOOLEAN isGB = FALSE;
2882
2883  /* compute a pertubed weight vector of the target weight vector */
2884  if(tp_deg > 1 && tp_deg <= nV) {
2885    ideal H0 = idHeadCC(G);
2886
2887    if (rParameter (currRing) != NULL)
2888      DefRingParlp();
2889    else
2890      VMrDefaultlp();
2891
2892    TargetRing = currRing;
2893    ssG = idrMoveR(G,EXXRing,currRing);
2894
2895    ideal H0_tmp = idrMoveR(H0,EXXRing,currRing);
2896    ideal H1 = idHeadCC(ssG);
2897
2898    /* Apply Lemma 2.2 in Collart et. al (1997) to check whether
2899       cone(k-1) is equal to cone(k) */
2900    if(test_G_GB_walk(H0_tmp,H1)==1) {
2901      idDelete(&H0_tmp);
2902      idDelete(&H1);
2903      G = ssG;
2904      ssG = NULL;
2905      newRing = currRing;
2906      delete ivNull;
2907
2908      if(npwinc != 0)
2909        goto LastGB_Finish;
2910      else {
2911        isGB = TRUE;
2912        goto KSTD_Finish;
2913      }
2914    }
2915    idDelete(&H0_tmp);
2916    idDelete(&H1);
2917
2918    iv_M_lp = MivMatrixOrderlp(nV);
2919    target_weight  = MPertVectors(ssG, iv_M_lp, tp_deg);
2920    delete iv_M_lp;
2921    //PrintS("\n// Input is not GB!!");
2922    rChangeCurrRing(EXXRing);
2923    G = idrMoveR(ssG, TargetRing,currRing);
2924
2925    if(Overflow_Error == TRUE)  {
2926      nOverflow_Error = Overflow_Error;
2927      NEG = 1;
2928      newRing = currRing;
2929      goto JUNI_STD;
2930    }
2931  }
2932
2933  while(1)
2934  {
2935    nwalk ++;
2936    nstep++;
2937
2938    if(nwalk==1)
2939      goto FIRST_STEP;
2940
2941    to=clock();
2942    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
2943    Gomega = MwalkInitialForm(G, curr_weight);
2944    xtif=xtif+clock()-to;
2945
2946#ifndef  BUCHBERGER_ALG
2947    if(isNolVector(curr_weight) == 0)
2948      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
2949    else
2950      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
2951#endif // BUCHBERGER_ALG
2952
2953    oldRing = currRing;
2954
2955    /* defiNe a new ring that its ordering is "(a(curr_weight),lp) */
2956    if (rParameter(currRing) != NULL)
2957      DefRingPar(curr_weight);
2958    else
2959      VMrDefault(curr_weight);
2960
2961    newRing = currRing;
2962    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
2963    to=clock();
2964    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
2965#ifdef  BUCHBERGER_ALG
2966    M = MstdhomCC(Gomega1);
2967#else
2968    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
2969    delete hilb_func;
2970#endif // BUCHBERGER_ALG
2971    xtstd=xtstd+clock()-to;
2972    /* change the ring to oldRing */
2973    rChangeCurrRing(oldRing);
2974    M1 =  idrMoveR(M, newRing,currRing);
2975    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
2976
2977     to=clock();
2978    /* compute a reduced Groebner basis of <G> w.r.t. "newRing" by the
2979       lifting process*/
2980    F = MLifttwoIdeal(Gomega2, M1, G);
2981    xtlift=xtlift+clock()-to;
2982    idDelete(&M1);
2983    idDelete(&Gomega2);
2984    idDelete(&G);
2985
2986    /* change the ring to newRing */
2987    rChangeCurrRing(newRing);
2988    F1 = idrMoveR(F, oldRing,currRing);
2989
2990    to=clock();
2991    /* reduce the Groebner basis <G> w.r.t. new ring */
2992    G = kInterRedCC(F1, NULL);
2993    xtred=xtred+clock()-to;
2994    idDelete(&F1);
2995
2996    if(endwalks == 1)
2997      break;
2998
2999  FIRST_STEP:
3000    to=clock();
3001    Overflow_Error = FALSE;
3002    /* compute a next weight vector */
3003    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
3004    xtnw=xtnw+clock()-to;
3005#ifdef PRINT_VECTORS
3006    MivString(curr_weight, target_weight, next_weight);
3007#endif
3008
3009    if(Overflow_Error == TRUE) {
3010      //PrintS("\n// ** The next vector does NOT stay in Cone!!\n");
3011#ifdef TEST_OVERFLOW
3012      goto  LastGB_Finish;
3013#endif
3014
3015      nnwinC = 0;
3016      if(tp_deg == nV)
3017        nlast = 1;
3018
3019      delete next_weight;
3020      break;
3021    }
3022
3023    if(MivComp(next_weight, ivNull) == 1) {
3024      //newRing = currRing;
3025      delete next_weight;
3026      break;
3027    }
3028
3029    if(MivComp(next_weight, target_weight) == 1)  {
3030      if(tp_deg == nV)
3031        endwalks = 1;
3032      else {
3033      REC_LAST_GB_ALT2:
3034        nOverflow_Error = Overflow_Error;
3035        tproc=tproc+clock()-tinput;
3036        /*
3037          Print("\n// takes %d steps and calls \"Rec_LastGB\" (%d):",
3038          nwalk, tp_deg+1);
3039        */
3040        G = Rec_LastGB(G,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
3041        newRing = currRing;
3042        delete next_weight;
3043        break;
3044      }
3045    }
3046
3047    for(i=nV-1; i>=0; i--) {
3048      //(*extra_curr_weight)[i] = (*curr_weight)[i];
3049      (*curr_weight)[i] = (*next_weight)[i];
3050    }
3051    delete next_weight;
3052  }//while
3053
3054  delete ivNull;
3055
3056  if(tp_deg != nV) {
3057    newRing = currRing;
3058
3059    if (rParameter(currRing) != NULL)
3060      DefRingParlp();
3061    else
3062      VMrDefaultlp();
3063
3064    F1 = idrMoveR(G, newRing,currRing);
3065
3066    if(nnwinC == 0 || test_w_in_ConeCC(F1, target_weight) != 1 ) {
3067      nOverflow_Error = Overflow_Error;
3068      //Print("\n//  takes %d steps and calls \"Rec_LastGB (%d):", tp_deg+1);
3069      tproc=tproc+clock()-tinput;
3070      F1 = Rec_LastGB(F1,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
3071    }
3072    delete target_weight;
3073
3074    TargetRing = currRing;
3075    rChangeCurrRing(EXXRing);
3076    result = idrMoveR(F1, TargetRing,currRing);
3077  }
3078  else  {
3079    if(nlast == 1) {
3080      JUNI_STD:
3081
3082      newRing = currRing;
3083      if (rParameter(currRing) != NULL)
3084        DefRingParlp();
3085      else
3086        VMrDefaultlp();
3087
3088      KSTD_Finish:
3089      if(isGB == FALSE)
3090        F1 = idrMoveR(G, newRing,currRing);
3091      else
3092        F1 = G;
3093      to=clock();
3094      // Print("\n// apply the Buchberger's alg in ring = %s",rString(currRing));
3095      // idElements(F1, "F1");
3096      G = MstdCC(F1);
3097      xtextra=xtextra+clock()-to;
3098
3099
3100      idDelete(&F1);
3101      newRing = currRing;
3102    }
3103
3104    LastGB_Finish:
3105    rChangeCurrRing(EXXRing);
3106    result = idrMoveR(G, newRing,currRing);
3107  }
3108
3109  if(Overflow_Error == FALSE)
3110    Overflow_Error=nError;
3111  /*
3112  Print("\n// \"Rec_LastGB\" (%d) took %d steps and %.2f sec.Overflow_Error (%d)", tp_deg,
3113        nwalk, ((double) tproc)/1000000, nOverflow_Error);
3114  */
3115  return(result);
3116}
3117
3118/* The following subroutine is the implementation of our second improved
3119   Groebner walk algorithm, i.e. the second altervative algorithm.
3120   First we use the Grobner walk algorithm and then we call the changed
3121   perturbation walk algorithm with increased degree, if an intermediate
3122   weight vector is equal to the current target weight vector.
3123   This call will be only repeated until we get the wanted reduced Groebner
3124   basis or n times, where n is the numbers of variables.
3125*/
3126/* walk + recursive LastGB */
3127ideal MAltwalk2(ideal Go, intvec* curr_weight, intvec* target_weight)
3128{
3129  Set_Error(FALSE);
3130  Overflow_Error = FALSE;
3131  BOOLEAN nOverflow_Error = FALSE;
3132  //Print("// pSetm_Error = (%d)", ErrorCheck());
3133
3134  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0;
3135  xftinput = clock();
3136  clock_t tostd, tproc;
3137
3138  nstep = 0;
3139  int i, nV = currRing->N;
3140  int nwalk=0, endwalks=0, nhilb=1;
3141
3142  ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G, G1;
3143  ring endRing, newRing, oldRing;
3144  intvec* ivNull = new intvec(nV);
3145  intvec* next_weight;
3146  intvec* extra_curr_weight = new intvec(nV);
3147  intvec* hilb_func;
3148  intvec* exivlp = Mivlp(nV);
3149
3150  ring XXRing = currRing;
3151
3152  //Print("\n// ring r_input = %s;", rString(currRing));
3153  to = clock();
3154  /* compute the reduced Groebner basis of the given ideal w.r.t.
3155     a "fast" monomial order, e.g. degree reverse lex. order (dp) */
3156  G = MstdCC(Go);
3157  tostd=clock()-to;
3158
3159  /*
3160  Print("\n// Computation of the first std took = %.2f sec",
3161        ((double) tostd)/1000000);
3162  */
3163  if(currRing->order[0] == ringorder_a)
3164    goto NEXT_VECTOR;
3165
3166  while(1)
3167  {
3168    nwalk ++;
3169    nstep ++;
3170    to = clock();
3171    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
3172    Gomega = MwalkInitialForm(G, curr_weight);
3173    xtif=xtif+clock()-to;
3174#if 0
3175    if(Overflow_Error == TRUE)
3176    {
3177      for(i=nV-1; i>=0; i--)
3178        (*curr_weight)[i] = (*extra_curr_weight)[i];
3179      delete extra_curr_weight;
3180      goto LAST_GB_ALT2;
3181    }
3182#endif
3183    oldRing = currRing;
3184
3185    /* define a new ring that its ordering is "(a(curr_weight),lp) */
3186    if (rParameter(currRing) != NULL)
3187      DefRingPar(curr_weight);
3188    else
3189      VMrDefault(curr_weight);
3190
3191    newRing = currRing;
3192    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
3193    to = clock();
3194    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
3195    M = MstdhomCC(Gomega1);
3196    xtstd=xtstd+clock()-to;
3197    /* change the ring to oldRing */
3198    rChangeCurrRing(oldRing);
3199    M1 =  idrMoveR(M, newRing,currRing);
3200    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
3201
3202    to = clock();
3203    /* compute the reduced Groebner basis of <G> w.r.t. "newRing"
3204       by the liftig process */
3205    F = MLifttwoIdeal(Gomega2, M1, G);
3206    xtlift=xtlift+clock()-to;
3207    idDelete(&M1);
3208    idDelete(&Gomega2);
3209    idDelete(&G);
3210
3211    /* change the ring to newRing */
3212    rChangeCurrRing(newRing);
3213    F1 = idrMoveR(F, oldRing,currRing);
3214
3215    to = clock();
3216    /* reduce the Groebner basis <G> w.r.t. newRing */
3217    G = kInterRedCC(F1, NULL);
3218    xtred=xtred+clock()-to;
3219    idDelete(&F1);
3220
3221    if(endwalks == 1)
3222      break;
3223
3224  NEXT_VECTOR:
3225    to = clock();
3226    /* compute a next weight vector */
3227    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
3228    xtnw=xtnw+clock()-to;
3229#ifdef PRINT_VECTORS
3230    MivString(curr_weight, target_weight, next_weight);
3231#endif
3232
3233    if(Overflow_Error == TRUE)
3234    {
3235      /*
3236        ivString(next_weight, "omega");
3237        PrintS("\n// ** The weight vector does NOT stay in Cone!!\n");
3238      */
3239#ifdef TEST_OVERFLOW
3240      goto  TEST_OVERFLOW_OI;
3241#endif
3242
3243
3244      newRing = currRing;
3245      if (rParameter(currRing) != NULL)
3246        DefRingPar(target_weight);
3247      else
3248        VMrDefault(target_weight);
3249
3250      F1 = idrMoveR(G, newRing,currRing);
3251      G = MstdCC(F1);
3252      idDelete(&F1);
3253      newRing = currRing;
3254      break;
3255    }
3256
3257    if(MivComp(next_weight, ivNull) == 1)
3258    {
3259      newRing = currRing;
3260      delete next_weight;
3261      break;
3262    }
3263
3264    if(MivComp(next_weight, target_weight) == 1)
3265    {
3266      if(MivSame(target_weight, exivlp)==1)
3267      {
3268      LAST_GB_ALT2:
3269        nOverflow_Error = Overflow_Error;
3270        tproc = clock()-xftinput;
3271        //Print("\n// takes %d steps and calls the recursion of level 2:",  nwalk);
3272        /* call the changed perturbation walk algorithm with degree 2 */
3273        G = Rec_LastGB(G, curr_weight, target_weight, 2,1);
3274        newRing = currRing;
3275        delete next_weight;
3276        break;
3277      }
3278      endwalks = 1;
3279    }
3280
3281    for(i=nV-1; i>=0; i--)
3282    {
3283      //(*extra_curr_weight)[i] = (*curr_weight)[i];
3284      (*curr_weight)[i] = (*next_weight)[i];
3285    }
3286    delete next_weight;
3287  }
3288 TEST_OVERFLOW_OI:
3289  rChangeCurrRing(XXRing);
3290  G = idrMoveR(G, newRing,currRing);
3291  delete ivNull;
3292  delete exivlp;
3293
3294#ifdef TIME_TEST
3295  Print("\n// \"Main procedure\"  took %d steps dnd %.2f sec. Overflow_Error (%d)", nwalk,
3296        ((double) tproc)/1000000, nOverflow_Error);
3297
3298  TimeStringFractal(xftinput, tostd, xtif, xtstd, xtextra,xtlift, xtred,xtnw);
3299
3300  Print("\n// pSetm_Error = (%d)", ErrorCheck());
3301  //Print("\n// Overflow_Error? (%d)", nOverflow_Error);
3302  Print("\n// Awalk2 took %d steps!!", nstep);
3303#endif
3304
3305  return(G);
3306}
3307
3308
3309/* 5.5.02 */
3310/* The implementation of the fractal walk algorithmus */
3311
3312/* The main procedur Mfwalk calls the recursive Subroutine
3313   rec_fractal_call to compute the wanted Grï¿œbner basis.
3314   At the main procedur we compute the reduced Grï¿œbner basis w.r.t. a "fast"
3315   order, e.g. "dp" and a sequence of weight vectors which are row vectors
3316   of a matrix. This matrix defines the given monomial order, e.g. "lp"
3317*/
3318
3319
3320
3321
3322/* perturb the matrix order of  "lex" */
3323static intvec* NewVectorlp(ideal I)
3324{
3325  int nV = currRing->N;
3326  intvec* iv_wlp =  MivMatrixOrderlp(nV);
3327  intvec* result = Mfpertvector(I, iv_wlp);
3328  delete iv_wlp;
3329  return result;
3330}
3331
3332int ngleich;
3333intvec* Xsigma;
3334intvec* Xtau;
3335int xn;
3336intvec* Xivinput;
3337intvec* Xivlp;
3338
3339
3340
3341/***********************************************************************
3342  The procedur REC_GB_Mwalk computes a GB for <G> w.r.t. the weight order
3343  otw, where G is a reduced GB w.r.t. the weight order cw.
3344  The new procedur Mwalk calls REC_GB.
3345************************************************************************/
3346static ideal REC_GB_Mwalk(ideal G, intvec* curr_weight, intvec* orig_target_weight,
3347                          int tp_deg, int npwinc)
3348{
3349  BOOLEAN nError = Overflow_Error;
3350  Overflow_Error = FALSE;
3351
3352  int i,  nV = currRing->N, ntwC,  npwinC;
3353  int nwalk=0, endwalks=0, nnwinC=1, nlast = 0;
3354  ideal Gomega, M, F, Gomega1, Gomega2, M1,F1,result,ssG;
3355  ring newRing, oldRing, TargetRing;
3356  intvec* iv_M_lp;
3357  intvec* target_weight;
3358  intvec* ivNull = new intvec(nV);
3359  intvec* hilb_func;
3360  BOOLEAN isGB = FALSE;
3361
3362  /* to avoid (1,0,...,0) as the target vector */
3363  intvec* last_omega = new intvec(nV);
3364  for(i=nV-1; i>0; i--)
3365    (*last_omega)[i] = 1;
3366  (*last_omega)[0] = 10000;
3367
3368  ring EXXRing = currRing;
3369
3370  /* compute a pertubed weight vector of the target weight vector */
3371  if(tp_deg > 1 && tp_deg <= nV)
3372  {
3373    ideal H0 = idHeadCC(G);
3374    if (rParameter(currRing) != NULL)
3375      DefRingPar(orig_target_weight);
3376    else
3377      VMrDefault(orig_target_weight);
3378
3379    TargetRing = currRing;
3380    ssG = idrMoveR(G,EXXRing,currRing);
3381
3382    ideal H0_tmp = idrMoveR(H0,EXXRing,currRing);
3383    ideal H1 = idHeadCC(ssG);
3384    id_Delete(&H0,EXXRing);
3385
3386    if(test_G_GB_walk(H0_tmp,H1)==1)
3387    {
3388      //Print("//input in %d-th recursive is a GB",tp_deg);
3389      idDelete(&H0_tmp);idDelete(&H1);
3390      G = ssG;
3391      ssG = NULL;
3392      newRing = currRing;
3393      delete ivNull;
3394      if(npwinc == 0)
3395      {
3396        isGB = TRUE;
3397        goto KSTD_Finish;
3398      }
3399      else
3400        goto LastGB_Finish;
3401    }
3402    idDelete(&H0_tmp);   idDelete(&H1);
3403
3404    intvec* ivlp = Mivlp(nV);
3405    if( MivSame(orig_target_weight, ivlp)==1 )
3406      iv_M_lp = MivMatrixOrderlp(nV);
3407    else
3408      iv_M_lp = MivMatrixOrder(orig_target_weight);
3409
3410    //target_weight  = MPertVectorslp(ssG, iv_M_lp, tp_deg);
3411    target_weight  = MPertVectors(ssG, iv_M_lp, tp_deg);
3412
3413    delete ivlp;
3414    delete iv_M_lp;
3415
3416    rChangeCurrRing(EXXRing);
3417    G = idrMoveR(ssG, TargetRing,currRing);
3418  }
3419
3420  while(1)
3421  {
3422    nwalk ++;
3423    nstep++;
3424    if(nwalk == 1)
3425      goto NEXT_STEP;
3426
3427    //Print("\n// Entering the %d-th step in the %d-th recursive:",nwalk,tp_deg);
3428    to = clock();
3429    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
3430    Gomega = MwalkInitialForm(G, curr_weight);
3431    xtif = xtif + clock()-to;
3432
3433#ifndef  BUCHBERGER_ALG
3434    if(isNolVector(curr_weight) == 0)
3435      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
3436    else
3437      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
3438#endif // BUCHBERGER_ALG
3439
3440    oldRing = currRing;
3441
3442    /* define a new ring that its ordering is "(a(curr_weight),lp) */
3443    if (rParameter(currRing) != NULL)
3444      DefRingPar(curr_weight);
3445    else
3446      VMrDefault(curr_weight);
3447
3448    newRing = currRing;
3449    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
3450
3451    to = clock();
3452    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
3453#ifdef  BUCHBERGER_ALG
3454    M = MstdhomCC(Gomega1);
3455#else
3456    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
3457    delete hilb_func;
3458#endif // BUCHBERGER_ALG
3459    xtstd = xtstd + clock() - to;
3460
3461    /* change the ring to oldRing */
3462    rChangeCurrRing(oldRing);
3463
3464    M1 =  idrMoveR(M, newRing,currRing);
3465    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
3466
3467    to = clock();
3468    F = MLifttwoIdeal(Gomega2, M1, G);
3469    xtlift = xtlift + clock() -to;
3470
3471    idDelete(&M1);
3472    idDelete(&Gomega2);
3473    idDelete(&G);
3474
3475
3476    /* change the ring to newRing */
3477    rChangeCurrRing(newRing);
3478    F1 = idrMoveR(F, oldRing,currRing);
3479
3480    to = clock();
3481    /* reduce the Groebner basis <G> w.r.t. new ring */
3482    G = kInterRedCC(F1, NULL);
3483    xtred = xtred + clock() -to;
3484
3485    idDelete(&F1);
3486
3487    if(endwalks == 1)
3488      break;
3489
3490  NEXT_STEP:
3491    to = clock();
3492    /* compute a next weight vector */
3493    intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
3494    xtnw = xtnw + clock() - to;
3495
3496#ifdef PRINT_VECTORS
3497    MivString(curr_weight, target_weight, next_weight);
3498#endif
3499
3500    /*check whether the computed vector does in the correct cone */
3501    //ntwC = test_w_in_ConeCC(G, next_weight);
3502    //if(ntwC != 1)
3503    if(Overflow_Error == TRUE)
3504    {
3505      PrintS("\n// ** The computed vector does NOT stay in Cone!!\n");
3506      nnwinC = 0;
3507      if(tp_deg == nV)
3508        nlast = 1;
3509      delete next_weight;
3510      break;
3511    }
3512    if(MivComp(next_weight, ivNull) == 1)
3513    {
3514      newRing = currRing;
3515      delete next_weight;
3516      break;
3517    }
3518
3519    if(MivComp(next_weight, target_weight) == 1)
3520    {
3521      if(tp_deg == nV)
3522        endwalks = 1;
3523      else {
3524        G = REC_GB_Mwalk(G,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
3525        newRing = currRing;
3526        delete next_weight;
3527        break;
3528      }
3529    }
3530
3531    /* 06.11.01 NOT Changed */
3532    for(i=nV-1; i>=0; i--)
3533      (*curr_weight)[i] = (*next_weight)[i];
3534
3535    delete next_weight;
3536  }//while
3537
3538  delete ivNull;
3539
3540  if(tp_deg != nV)
3541  {
3542    //28.07.03
3543    newRing = currRing;
3544
3545    if (rParameter(currRing) != NULL)
3546      //  DefRingParlp(); //
3547      DefRingPar(orig_target_weight);
3548    else
3549      VMrDefault(orig_target_weight);
3550
3551
3552    F1 = idrMoveR(G, newRing,currRing);
3553
3554    if(nnwinC == 0)
3555      F1 = REC_GB_Mwalk(F1,curr_weight, orig_target_weight, tp_deg+1,nnwinC);
3556    else
3557      if(test_w_in_ConeCC(F1, target_weight) != 1)
3558        F1 = REC_GB_Mwalk(F1,curr_weight, orig_target_weight,tp_deg+1,nnwinC);
3559
3560    delete target_weight;
3561
3562    TargetRing = currRing;
3563    rChangeCurrRing(EXXRing);
3564    result = idrMoveR(F1, TargetRing,currRing);
3565  }
3566  else
3567  {
3568    if(nlast == 1)
3569    {
3570      if (rParameter(currRing) != NULL)
3571        DefRingPar(orig_target_weight);
3572      else
3573        VMrDefault(orig_target_weight);
3574
3575
3576    KSTD_Finish:
3577      if(isGB == FALSE)
3578        F1 = idrMoveR(G, newRing,currRing);
3579      else
3580        F1 = G;
3581      to=clock();
3582      /* apply Buchberger alg to compute a red. GB of F1 */
3583      G = MstdCC(F1);
3584      xtextra=clock()-to;
3585      idDelete(&F1);
3586      newRing = currRing;
3587    }
3588
3589  LastGB_Finish:
3590    rChangeCurrRing(EXXRing);
3591    result = idrMoveR(G, newRing,currRing);
3592  }
3593
3594  if(Overflow_Error == FALSE)
3595    Overflow_Error = nError;
3596
3597  return(result);
3598}
3599
3600
3601/* 08.09.02 */
3602/******** THE NEW GRï¿œBNER WALK ALGORITHM **********/
3603/* Grï¿œbnerwalk with a recursive "second" alternative GW, REC_GB_Mwalk
3604   that only computes the last reduced GB */
3605ideal Mwalk(ideal Go, intvec* curr_weight, intvec* target_weight)
3606{
3607  Set_Error(FALSE);
3608  Overflow_Error = FALSE;
3609  //Print("// pSetm_Error = (%d)", ErrorCheck());
3610
3611  clock_t tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
3612  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
3613  tinput = clock();
3614  clock_t tim;
3615  nstep=0;
3616  int i, nV = currRing->N;
3617  int nwalk=0, endwalks=0;
3618
3619  ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G, G1;
3620  ring endRing, newRing, oldRing;
3621  intvec* ivNull = new intvec(nV);
3622  intvec* exivlp = Mivlp(nV);
3623  intvec* hilb_func;
3624
3625  intvec* tmp_weight = new intvec(nV);
3626  for(i=nV-1; i>=0; i--)
3627    (*tmp_weight)[i] = (*curr_weight)[i];
3628
3629   /* to avoid (1,0,...,0) as the target vector */
3630  intvec* last_omega = new intvec(nV);
3631  for(i=nV-1; i>0; i--)
3632    (*last_omega)[i] = 1;
3633  (*last_omega)[0] = 10000;
3634
3635  ring XXRing = currRing;
3636
3637  to = clock();
3638  /* the monomial ordering of this current ring would be "dp" */
3639  G = MstdCC(Go);
3640  tostd = clock()-to;
3641
3642  if(currRing->order[0] == ringorder_a)
3643    goto NEXT_VECTOR;
3644
3645  while(1)
3646  {
3647    nwalk ++;
3648    nstep ++;
3649    to = clock();
3650    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
3651    Gomega = MwalkInitialForm(G, curr_weight);
3652    tif = tif + clock()-to;
3653    oldRing = currRing;
3654
3655    if(endwalks == 1)
3656    {
3657      /* compute a reduced Groebner basis of Gomega w.r.t. >>_cw by
3658         the recursive changed perturbation walk alg. */
3659      tim = clock();
3660      /*
3661        Print("\n// **** Grï¿œbnerwalk took %d steps and ", nwalk);
3662        PrintS("\n// **** call the rec. Pert. Walk to compute a red GB of:");
3663        idElements(Gomega, "G_omega");
3664      */
3665
3666      if(MivSame(exivlp, target_weight)==1)
3667        M = REC_GB_Mwalk(idCopy(Gomega), tmp_weight, curr_weight, 2,1);
3668      else
3669        goto NORMAL_GW;
3670      /*
3671        Print("\n//  time for the last std(Gw)  = %.2f sec",
3672        ((double) (clock()-tim)/1000000));
3673        PrintS("\n// ***************************************************\n");
3674      */
3675#ifdef CHECK_IDEAL_MWALK
3676      idElements(Gomega, "G_omega");
3677      headidString(Gomega, "Gw");
3678      idElements(M, "M");
3679      //headidString(M, "M");
3680#endif
3681      to = clock();
3682      F = MLifttwoIdeal(Gomega, M, G);
3683      xtlift = xtlift + clock() - to;
3684
3685      idDelete(&Gomega);
3686      idDelete(&M);
3687      idDelete(&G);
3688
3689      oldRing = currRing;
3690
3691      /* create a new ring newRing */
3692       if (rParameter(currRing) != NULL)
3693         DefRingPar(curr_weight);
3694       else
3695         VMrDefault(curr_weight);
3696
3697      newRing = currRing;
3698      F1 = idrMoveR(F, oldRing,currRing);
3699    }
3700    else
3701    {
3702    NORMAL_GW:
3703#ifndef  BUCHBERGER_ALG
3704      if(isNolVector(curr_weight) == 0)
3705        hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
3706      else
3707        hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
3708#endif // BUCHBERGER_ALG
3709
3710      /* define a new ring that its ordering is "(a(curr_weight),lp) */
3711      if (rParameter(currRing) != NULL)
3712        DefRingPar(curr_weight);
3713      else
3714        VMrDefault(curr_weight);
3715
3716      newRing = currRing;
3717      Gomega1 = idrMoveR(Gomega, oldRing,currRing);
3718
3719      to = clock();
3720      /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
3721#ifdef  BUCHBERGER_ALG
3722      M = MstdhomCC(Gomega1);
3723#else
3724      M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
3725      delete hilb_func;
3726#endif // BUCHBERGER_ALG
3727      tstd = tstd + clock() - to;
3728
3729      /* change the ring to oldRing */
3730      rChangeCurrRing(oldRing);
3731      M1 =  idrMoveR(M, newRing,currRing);
3732      Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
3733
3734      to = clock();
3735      /* compute a representation of the generators of submod (M)
3736         with respect to those of mod (Gomega).
3737         Gomega is a reduced Groebner basis w.r.t. the current ring */
3738      F = MLifttwoIdeal(Gomega2, M1, G);
3739      tlift = tlift + clock() - to;
3740
3741      idDelete(&M1);
3742      idDelete(&Gomega2);
3743      idDelete(&G);
3744
3745      /* change the ring to newRing */
3746      rChangeCurrRing(newRing);
3747      F1 = idrMoveR(F, oldRing,currRing);
3748    }
3749
3750    to = clock();
3751    /* reduce the Groebner basis <G> w.r.t. new ring */
3752    G = kInterRedCC(F1, NULL);
3753    if(endwalks != 1)
3754      tred = tred + clock() - to;
3755    else
3756      xtred = xtred + clock() - to;
3757
3758    idDelete(&F1);
3759    if(endwalks == 1)
3760      break;
3761
3762  NEXT_VECTOR:
3763    to = clock();
3764    /* compute a next weight vector */
3765    intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G);
3766    tnw = tnw + clock() - to;
3767#ifdef PRINT_VECTORS
3768    MivString(curr_weight, target_weight, next_weight);
3769#endif
3770
3771    //if(test_w_in_ConeCC(G, next_weight) != 1)
3772    if(Overflow_Error == TRUE)
3773    {
3774      newRing = currRing;
3775      PrintS("\n// ** The computed vector does NOT stay in Cone!!\n");
3776
3777      if (rParameter(currRing) != NULL)
3778        DefRingPar(target_weight);
3779      else
3780        VMrDefault(target_weight);
3781
3782      F1 = idrMoveR(G, newRing,currRing);
3783      G = MstdCC(F1);
3784      idDelete(&F1);
3785
3786      newRing = currRing;
3787      break;
3788    }
3789
3790    if(MivComp(next_weight, ivNull) == 1)
3791    {
3792      newRing = currRing;
3793      delete next_weight;
3794      break;
3795    }
3796    if(MivComp(next_weight, target_weight) == 1)
3797      endwalks = 1;
3798
3799    for(i=nV-1; i>=0; i--) {
3800      (*tmp_weight)[i] = (*curr_weight)[i];
3801      (*curr_weight)[i] = (*next_weight)[i];
3802    }
3803    delete next_weight;
3804  }
3805  rChangeCurrRing(XXRing);
3806  G = idrMoveR(G, newRing,currRing);
3807
3808  delete tmp_weight;
3809  delete ivNull;
3810  delete exivlp;
3811
3812#ifdef TIME_TEST
3813  TimeString(tinput, tostd, tif, tstd, tlift, tred, tnw, nstep);
3814
3815  Print("\n// pSetm_Error = (%d)", ErrorCheck());
3816  Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
3817#endif
3818  return(G);
3819}
3820
3821  /* 2.12.02*/
3822ideal Mwalk_tst(ideal Go, intvec* curr_weight, intvec* target_weight)
3823{
3824  clock_t tinput=clock();
3825  //idString(Go,"Ginp");
3826  int i, nV = currRing->N;
3827  int nwalk=0, endwalks=0;
3828
3829  ideal Gomega, M, F, Gomega1, Gomega2, M1, F1, G, G1;
3830  ring endRing, newRing, oldRing;
3831  intvec* ivNull = new intvec(nV);
3832  ring XXRing = currRing;
3833
3834  intvec* tmp_weight = new intvec(nV);
3835  for(i=nV-1; i>=0; i--)
3836    (*tmp_weight)[i] = (*curr_weight)[i];
3837
3838  /* the monomial ordering of this current ring would be "dp" */
3839  G = MstdCC(Go);
3840
3841  intvec* hilb_func;
3842
3843  /* to avoid (1,0,...,0) as the target vector */
3844  intvec* last_omega = new intvec(nV);
3845  for(i=nV-1; i>0; i--)
3846    (*last_omega)[i] = 1;
3847  (*last_omega)[0] = 10000;
3848
3849  while(1)
3850  {
3851    nwalk ++;
3852    //Print("\n// Entering the %d-th step:", nwalk);
3853    //Print("\n// ring r[%d] = %s;", nwalk, rString(currRing));
3854    idString(G,"G");
3855    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
3856    Gomega = MwalkInitialForm(G, curr_weight);
3857    //ivString(curr_weight, "omega");
3858    idString(Gomega,"Gw");
3859
3860#ifndef  BUCHBERGER_ALG
3861    if(isNolVector(curr_weight) == 0)
3862      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
3863    else
3864      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
3865#endif // BUCHBERGER_ALG
3866
3867
3868    oldRing = currRing;
3869
3870    /* define a new ring that its ordering is "(a(curr_weight),lp) */
3871    VMrDefault(curr_weight);
3872    newRing = currRing;
3873
3874    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
3875
3876    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
3877#ifdef  BUCHBERGER_ALG
3878    M = MstdhomCC(Gomega1);
3879#else
3880    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
3881    delete hilb_func;
3882#endif // BUCHBERGER_ALG
3883
3884    idString(M,"M");
3885
3886      /* change the ring to oldRing */
3887    rChangeCurrRing(oldRing);
3888    M1 =  idrMoveR(M, newRing,currRing);
3889    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
3890
3891      /* compute a representation of the generators of submod (M)
3892         with respect to those of mod (Gomega).
3893         Gomega is a reduced Groebner basis w.r.t. the current ring */
3894    F = MLifttwoIdeal(Gomega2, M1, G);
3895    idDelete(&M1);
3896    idDelete(&Gomega2);
3897    idDelete(&G);
3898    idString(F,"F");
3899
3900    /* change the ring to newRing */
3901    rChangeCurrRing(newRing);
3902    F1 = idrMoveR(F, oldRing,currRing);
3903
3904    /* reduce the Groebner basis <G> w.r.t. new ring */
3905    G = kInterRedCC(F1, NULL);
3906    //idSkipZeroes(G);//done by kInterRed
3907    idDelete(&F1);
3908    idString(G,"G");
3909    if(endwalks == 1)
3910      break;
3911
3912    /* compute a next weight vector */
3913    intvec* next_weight = MkInterRedNextWeight(curr_weight,target_weight,G);
3914#ifdef PRINT_VECTORS
3915    MivString(curr_weight, target_weight, next_weight);
3916#endif
3917
3918    if(MivComp(next_weight, ivNull) == 1)
3919    {
3920      delete next_weight;
3921      break;
3922    }
3923    if(MivComp(next_weight, target_weight) == 1)
3924      endwalks = 1;
3925
3926    for(i=nV-1; i>=0; i--)
3927      (*tmp_weight)[i] = (*curr_weight)[i];
3928
3929    /* 06.11.01 to free the memory: NOT Changed!!*/
3930    for(i=nV-1; i>=0; i--)
3931      (*curr_weight)[i] = (*next_weight)[i];
3932    delete next_weight;
3933  }
3934  rChangeCurrRing(XXRing);
3935  G = idrMoveR(G, newRing,currRing);
3936
3937  delete tmp_weight;
3938  delete ivNull;
3939  PrintLn();
3940  return(G);
3941}
3942
3943
3944
3945/**************************************************************/
3946/*     Implementation of the perturbation walk algorithm      */
3947/**************************************************************/
3948/* If the perturbed target weight vector or an intermediate weight vector
3949   doesn't stay in the correct Groebner cone, we have only
3950   a reduced Groebner basis for the given ideal with respect to
3951   a monomial order which differs to the given order.
3952   Then we have to compute the wanted  reduced Groebner basis for it.
3953   For this, we can use
3954   1) the improved Buchberger algorithm or
3955   2) the changed perturbation walk algorithm with a decreased degree.
3956*/
3957/* use kStd, if nP = 0, else call LastGB */
3958ideal Mpwalk(ideal Go, int op_deg, int tp_deg,intvec* curr_weight,
3959             intvec* target_weight, int nP)
3960{
3961  Set_Error(FALSE  );
3962  Overflow_Error = FALSE;
3963  //Print("// pSetm_Error = (%d)", ErrorCheck());
3964
3965  clock_t  tinput, tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0;
3966  xtextra=0;
3967  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0;
3968  tinput = clock();
3969
3970  clock_t tim;
3971
3972  nstep = 0;
3973  int i, ntwC=1, ntestw=1,  nV = currRing->N, op_tmp = op_deg;
3974  int endwalks=0, nhilb=0, ntestomega=0;
3975
3976  ideal Gomega, M, F, G, Gomega1, Gomega2, M1,F1,Eresult,ssG;
3977  ring newRing, oldRing, TargetRing;
3978  intvec* iv_M_dp;
3979  intvec* iv_M_lp;
3980  intvec* exivlp = Mivlp(nV);
3981  intvec* orig_target = target_weight;
3982  intvec* pert_target_vector = target_weight;
3983  intvec* ivNull = new intvec(nV);
3984  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
3985  intvec* cw_tmp = curr_weight;
3986  intvec* hilb_func;
3987  intvec* next_weight;
3988  intvec* extra_curr_weight = new intvec(nV);
3989
3990  /* to avoid (1,0,...,0) as the target vector */
3991  intvec* last_omega = new intvec(nV);
3992  for(i=nV-1; i>0; i--)
3993    (*last_omega)[i] = 1;
3994  (*last_omega)[0] = 10000;
3995
3996  ring XXRing = currRing;
3997
3998
3999  to = clock();
4000  /* perturbs the original vector */
4001  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) := "dp"
4002  {
4003    G = MstdCC(Go);
4004    tostd = clock()-to;
4005    if(op_deg != 1){
4006      iv_M_dp = MivMatrixOrderdp(nV);
4007      //ivString(iv_M_dp, "iv_M_dp");
4008      curr_weight = MPertVectors(G, iv_M_dp, op_deg);
4009    }
4010  }
4011  else
4012  {
4013    //ring order := (a(curr_weight),lp);
4014    if (rParameter(currRing) != NULL)
4015      DefRingPar(curr_weight);
4016    else
4017      VMrDefault(curr_weight);
4018
4019    G = idrMoveR(Go, XXRing,currRing);
4020    G = MstdCC(G);
4021    tostd = clock()-to;
4022    if(op_deg != 1){
4023      iv_M_dp = MivMatrixOrder(curr_weight);
4024      curr_weight = MPertVectors(G, iv_M_dp, op_deg);
4025    }
4026  }
4027  delete iv_dp;
4028  if(op_deg != 1) delete iv_M_dp;
4029
4030  ring HelpRing = currRing;
4031
4032  /* perturbs the target weight vector */
4033  if(tp_deg > 1 && tp_deg <= nV)
4034  {
4035    if (rParameter(currRing) != NULL)
4036      DefRingPar(target_weight);
4037    else
4038      VMrDefault(target_weight);
4039
4040    TargetRing = currRing;
4041    ssG = idrMoveR(G,HelpRing,currRing);
4042    if(MivSame(target_weight, exivlp) == 1)
4043    {
4044      iv_M_lp = MivMatrixOrderlp(nV);
4045      //ivString(iv_M_lp, "iv_M_lp");
4046      //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
4047      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
4048    }
4049    else
4050    {
4051      iv_M_lp = MivMatrixOrder(target_weight);
4052      //target_weight = MPertVectorslp(ssG, iv_M_lp, tp_deg);
4053      target_weight = MPertVectors(ssG, iv_M_lp, tp_deg);
4054    }
4055    delete iv_M_lp;
4056    pert_target_vector = target_weight; //vor 19. mai 2003//test 19 Junu 03
4057    rChangeCurrRing(HelpRing);
4058    G = idrMoveR(ssG, TargetRing,currRing);
4059  }
4060  /*
4061    Print("\n// Perturbationwalkalg. vom Gradpaar (%d,%d):",op_deg,tp_deg);
4062    ivString(curr_weight, "new sigma");
4063    ivString(target_weight, "new tau");
4064  */
4065  while(1)
4066  {
4067    nstep ++;
4068    to = clock();
4069    /* compute an initial form ideal of <G> w.r.t. the weight vector
4070       "curr_weight" */
4071    Gomega = MwalkInitialForm(G, curr_weight);
4072
4073
4074#ifdef ENDWALKS
4075    if(endwalks == 1){
4076      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
4077      idElements(G, "G");
4078      // idElements(Gomega, "Gw");
4079      headidString(G, "G");
4080      //headidString(Gomega, "Gw");
4081    }
4082#endif
4083
4084    tif = tif + clock()-to;
4085
4086#ifndef  BUCHBERGER_ALG
4087    if(isNolVector(curr_weight) == 0)
4088      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
4089    else
4090      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
4091#endif // BUCHBERGER_ALG
4092
4093    oldRing = currRing;
4094
4095    /* define a new ring that its ordering is "(a(curr_weight),lp) */
4096    if (rParameter(currRing) != NULL)
4097      DefRingPar(curr_weight);
4098    else
4099      VMrDefault(curr_weight);
4100
4101    newRing = currRing;
4102    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
4103
4104#ifdef ENDWALKS
4105    if(endwalks==1)
4106    {
4107      Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
4108      idElements(Gomega1, "Gw");
4109      headidString(Gomega1, "headGw");
4110      PrintS("\n// compute a rGB of Gw:\n");
4111
4112#ifndef  BUCHBERGER_ALG
4113      ivString(hilb_func, "w");
4114#endif
4115    }
4116#endif
4117
4118    tim = clock();
4119    to = clock();
4120    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
4121#ifdef  BUCHBERGER_ALG
4122    M = MstdhomCC(Gomega1);
4123#else
4124    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
4125    delete hilb_func;
4126#endif // BUCHBERGER_ALG
4127
4128    if(endwalks == 1){
4129      xtstd = xtstd+clock()-to;
4130#ifdef ENDWALKS
4131      Print("\n// time for the last std(Gw)  = %.2f sec\n",
4132            ((double) clock())/1000000 -((double)tim) /1000000);
4133#endif
4134    }
4135    else
4136      tstd=tstd+clock()-to;
4137
4138    /* change the ring to oldRing */
4139    rChangeCurrRing(oldRing);
4140    M1 =  idrMoveR(M, newRing,currRing);
4141    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
4142
4143    //if(endwalks==1)  PrintS("\n// Lifting is working:..");
4144
4145    to=clock();
4146    /* compute a representation of the generators of submod (M)
4147       with respect to those of mod (Gomega).
4148       Gomega is a reduced Groebner basis w.r.t. the current ring */
4149    F = MLifttwoIdeal(Gomega2, M1, G);
4150    if(endwalks != 1)
4151      tlift = tlift+clock()-to;
4152    else
4153      xtlift=clock()-to;
4154
4155    idDelete(&M1);
4156    idDelete(&Gomega2);
4157    idDelete(&G);
4158
4159    /* change the ring to newRing */
4160    rChangeCurrRing(newRing);
4161    F1 = idrMoveR(F, oldRing,currRing);
4162
4163    //if(endwalks==1)PrintS("\n// InterRed is working now:");
4164
4165    to=clock();
4166    /* reduce the Groebner basis <G> w.r.t. new ring */
4167    G = kInterRedCC(F1, NULL);
4168    if(endwalks != 1)
4169      tred = tred+clock()-to;
4170    else
4171      xtred=clock()-to;
4172
4173    idDelete(&F1);
4174
4175    if(endwalks == 1)
4176      break;
4177
4178    to=clock();
4179    /* compute a next weight vector */
4180    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
4181    tnw=tnw+clock()-to;
4182#ifdef PRINT_VECTORS
4183    MivString(curr_weight, target_weight, next_weight);
4184#endif
4185
4186    if(Overflow_Error == TRUE)
4187    {
4188      ntwC = 0;
4189      ntestomega = 1;
4190      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
4191      //idElements(G, "G");
4192      delete next_weight;
4193      goto FINISH_160302;
4194    }
4195    if(MivComp(next_weight, ivNull) == 1){
4196      newRing = currRing;
4197      delete next_weight;//16.03.02
4198      //Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
4199      break;
4200    }
4201    if(MivComp(next_weight, target_weight) == 1)
4202      endwalks = 1;
4203
4204    for(i=nV-1; i>=0; i--)
4205      (*curr_weight)[i] = (*next_weight)[i];
4206
4207    delete next_weight;
4208  }//while
4209
4210  if(tp_deg != 1)
4211  {
4212  FINISH_160302://16.03.02
4213    if(MivSame(orig_target, exivlp) == 1)
4214      if (rParameter(currRing) != NULL)
4215        DefRingParlp();
4216      else
4217        VMrDefaultlp();
4218    else
4219      if (rParameter(currRing) != NULL)
4220        DefRingPar(orig_target);
4221      else
4222        VMrDefault(orig_target);
4223
4224    TargetRing=currRing;
4225    F1 = idrMoveR(G, newRing,currRing);
4226#ifdef CHECK_IDEAL
4227      headidString(G, "G");
4228#endif
4229
4230
4231    // check whether the pertubed target vector stays in the correct cone
4232    if(ntwC != 0){
4233      ntestw = test_w_in_ConeCC(F1, pert_target_vector);
4234    }
4235
4236    if( ntestw != 1 || ntwC == 0)
4237    {
4238      /*
4239      if(ntestw != 1){
4240        ivString(pert_target_vector, "tau");
4241        PrintS("\n// ** perturbed target vector doesn't stay in cone!!");
4242        Print("\n// ring r%d = %s;\n", nstep, rString(currRing));
4243        idElements(F1, "G");
4244      }
4245      */
4246      /* LastGB is "better" than the kStd subroutine */
4247      to=clock();
4248      ideal eF1;
4249      if(nP == 0 || tp_deg == 1 || MivSame(orig_target, exivlp) != 1){
4250        // PrintS("\n// ** calls \"std\" to compute a GB");
4251        eF1 = MstdCC(F1);
4252        idDelete(&F1);
4253      }
4254      else {
4255        // PrintS("\n// ** calls \"LastGB\" to compute a GB");
4256        rChangeCurrRing(newRing);
4257        ideal F2 = idrMoveR(F1, TargetRing,currRing);
4258        eF1 = LastGB(F2, curr_weight, tp_deg-1);
4259        F2=NULL;
4260      }
4261      xtextra=clock()-to;
4262      ring exTargetRing = currRing;
4263
4264      rChangeCurrRing(XXRing);
4265      Eresult = idrMoveR(eF1, exTargetRing,currRing);
4266    }
4267    else{
4268      rChangeCurrRing(XXRing);
4269      Eresult = idrMoveR(F1, TargetRing,currRing);
4270    }
4271  }
4272  else {
4273    rChangeCurrRing(XXRing);
4274    Eresult = idrMoveR(G, newRing,currRing);
4275  }
4276  delete ivNull;
4277  if(tp_deg != 1)
4278    delete target_weight;
4279
4280  if(op_deg != 1 )
4281    delete curr_weight;
4282
4283  delete exivlp;
4284  delete last_omega;
4285
4286#ifdef TIME_TEST
4287  TimeStringFractal(tinput, tostd, tif+xtif, tstd+xtstd,0, tlift+xtlift, tred+xtred,
4288             tnw+xtnw);
4289
4290  Print("\n// pSetm_Error = (%d)", ErrorCheck());
4291  Print("\n// It took %d steps and Overflow_Error? (%d)\n", nstep,  Overflow_Error);
4292#endif
4293  return(Eresult);
4294}
4295
4296intvec* XivNull;
4297
4298/* define a matrix (1 ... 1) */
4299intvec* MMatrixone(int nV)
4300{
4301  int i,j;
4302  intvec* ivM = new intvec(nV*nV);
4303
4304  for(i=0; i<nV; i++)
4305    for(j=0; j<nV; j++)
4306    (*ivM)[i*nV + j] = 1;
4307
4308  return(ivM);
4309}
4310
4311int nnflow;
4312int Xcall;
4313int Xngleich;
4314
4315/* 27.07.03 */
4316/* Perturb the start weight vector at the top level, i.e. nlev = 1 */
4317static ideal rec_fractal_call(ideal G, int nlev, intvec* omtmp)
4318{
4319  Overflow_Error =  FALSE;
4320  //Print("\n\n// Entering the %d-th recursion:", nlev);
4321
4322  int i, nV = currRing->N;
4323  ring new_ring, testring, extoRing;
4324  ideal Gomega, Gomega1, Gomega2, F, F1, Gresult, Gresult1, G1, Gt;
4325  int nwalks = 0;
4326  intvec* Mwlp;
4327  intvec* hilb_func;
4328  intvec* extXtau;
4329  intvec* next_vect;
4330  intvec* omega2 = new intvec(nV);
4331  intvec* altomega = new intvec(nV);
4332
4333  BOOLEAN isnewtarget = FALSE;
4334
4335  /* to avoid (1,0,...,0) as the target vector (Hans) */
4336  intvec* last_omega = new intvec(nV);
4337  for(i=nV-1; i>0; i--)
4338    (*last_omega)[i] = 1;
4339  (*last_omega)[0] = 10000;
4340
4341  intvec* omega = new intvec(nV);
4342  for(i=0; i<nV; i++) {
4343    if(Xsigma->length() == nV)
4344      (*omega)[i] =  (*Xsigma)[i];
4345    else
4346      (*omega)[i] = (*Xsigma)[(nV*(nlev-1))+i];
4347
4348    (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
4349  }
4350
4351   if(nlev == 1)  Xcall = 1;
4352   else Xcall = 0;
4353
4354  ring oRing = currRing;
4355
4356  while(1)
4357  {
4358#ifdef FIRST_STEP_FRACTAL
4359    // perturb the current weight vector only on the top level or
4360    // after perturbation of the both vectors, nlev = 2 as the top level
4361    if((nlev == 1 && Xcall == 0) || (nlev == 2 && Xngleich == 1))
4362      if(islengthpoly2(G) == 1)
4363      {
4364        Mwlp = MivWeightOrderlp(omega);
4365        Xsigma = Mfpertvector(G, Mwlp);
4366        delete Mwlp;
4367        Overflow_Error = FALSE;
4368      }
4369#endif
4370    nwalks ++;
4371    NEXT_VECTOR_FRACTAL:
4372    to=clock();
4373    /* determine the next border */
4374    next_vect = MkInterRedNextWeight(omega,omega2,G);
4375    xtnw=xtnw+clock()-to;
4376#ifdef PRINT_VECTORS
4377    MivString(omega, omega2, next_vect);
4378#endif
4379    oRing = currRing;
4380
4381    /* We only perturb the current target vector at the recursion  level 1 */
4382    if(Xngleich == 0 && nlev == 1) //(ngleich == 0) important, e.g. ex2, ex3
4383      if (MivComp(next_vect, omega2) == 1)
4384      {
4385        /* to dispense with taking initial (and lifting/interreducing
4386           after the call of recursion */
4387        //Print("\n\n// ** Perturb the both vectors with degree %d with",nlev);
4388        //idElements(G, "G");
4389
4390        Xngleich = 1;
4391        nlev +=1;
4392
4393        if (rParameter(currRing) != NULL)
4394          DefRingPar(omtmp);
4395        else
4396          VMrDefault(omtmp);
4397
4398        testring = currRing;
4399        Gt = idrMoveR(G, oRing,currRing);
4400
4401        /* perturb the original target vector w.r.t. the current GB */
4402        delete Xtau;
4403        Xtau = NewVectorlp(Gt);
4404
4405        rChangeCurrRing(oRing);
4406        G = idrMoveR(Gt, testring,currRing);
4407
4408        /* perturb the current vector w.r.t. the current GB */
4409        Mwlp = MivWeightOrderlp(omega);
4410        Xsigma = Mfpertvector(G, Mwlp);
4411        delete Mwlp;
4412
4413        for(i=nV-1; i>=0; i--) {
4414          (*omega2)[i] = (*Xtau)[nV+i];
4415          (*omega)[i] = (*Xsigma)[nV+i];
4416        }
4417
4418        delete next_vect;
4419        to=clock();
4420
4421        /* to avoid the value of Overflow_Error that occur in Mfpertvector*/
4422        Overflow_Error = FALSE;
4423
4424        next_vect = MkInterRedNextWeight(omega,omega2,G);
4425        xtnw=xtnw+clock()-to;
4426
4427#ifdef PRINT_VECTORS
4428        MivString(omega, omega2, next_vect);
4429#endif
4430      }
4431
4432
4433    /* check whether the the computed vector is in the correct cone */
4434    /* If no, the reduced GB of an omega-homogeneous ideal will be
4435       computed by Buchberger algorithm and stop this recursion step*/
4436    //if(test_w_in_ConeCC(G, next_vect) != 1) //e.g. Example s7, cyc6
4437    if(Overflow_Error == TRUE)
4438    {
4439      delete next_vect;
4440
4441    OVERFLOW_IN_NEXT_VECTOR:
4442      if (rParameter(currRing) != NULL)
4443        DefRingPar(omtmp);
4444      else
4445        VMrDefault(omtmp);
4446
4447#ifdef TEST_OVERFLOW
4448      Gt = idrMoveR(G, oRing,currRing);
4449      Gt = NULL; return(Gt);
4450#endif
4451
4452      //Print("\n\n// apply BB's alg. in ring r = %s;", rString(currRing));
4453      to=clock();
4454      Gt = idrMoveR(G, oRing,currRing);
4455      G1 = MstdCC(Gt);
4456      xtextra=xtextra+clock()-to;
4457      Gt = NULL;
4458
4459      delete omega2;
4460      delete altomega;
4461
4462      //Print("\n// Leaving the %d-th recursion with %d steps", nlev, nwalks);
4463      //Print("  ** Overflow_Error? (%d)", Overflow_Error);
4464      nnflow ++;
4465
4466      Overflow_Error = FALSE;
4467      return (G1);
4468    }
4469
4470
4471    /* If the perturbed target vector stays in the correct cone,
4472       return the current GB,
4473       otherwise, return the computed  GB by the Buchberger-algorithm.
4474       Then we update the perturbed target vectors w.r.t. this GB. */
4475
4476    /* the computed vector is equal to the origin vector, since
4477       t is not defined */
4478    if (MivComp(next_vect, XivNull) == 1)
4479    {
4480      if (rParameter(currRing) != NULL)
4481        DefRingPar(omtmp);
4482      else
4483        VMrDefault(omtmp);
4484
4485      testring = currRing;
4486      Gt = idrMoveR(G, oRing,currRing);
4487
4488      if(test_w_in_ConeCC(Gt, omega2) == 1) {
4489        delete omega2;
4490        delete next_vect;
4491        delete altomega;
4492        //Print("\n// Leaving the %d-th recursion with %d steps ",nlev, nwalks);
4493        //Print(" ** Overflow_Error? (%d)", Overflow_Error);
4494
4495        return (Gt);
4496      }
4497      else
4498      {
4499        //ivString(omega2, "tau'");
4500        //Print("\n//  tau' doesn't stay in the correct cone!!");
4501
4502#ifndef  MSTDCC_FRACTAL
4503        //07.08.03
4504        //ivString(Xtau, "old Xtau");
4505        intvec* Xtautmp = Mfpertvector(Gt, MivMatrixOrder(omtmp));
4506#ifdef TEST_OVERFLOW
4507      if(Overflow_Error == TRUE)
4508      Gt = NULL; return(Gt);
4509#endif
4510
4511        if(MivSame(Xtau, Xtautmp) == 1)
4512        {
4513          //PrintS("\n// Update vectors are equal to the old vectors!!");
4514          delete Xtautmp;
4515          goto FRACTAL_MSTDCC;
4516        }
4517
4518        Xtau = Xtautmp;
4519        Xtautmp = NULL;
4520        //ivString(Xtau, "new  Xtau");
4521
4522        for(i=nV-1; i>=0; i--)
4523          (*omega2)[i] = (*Xtau)[(nlev-1)*nV+i];
4524
4525        //Print("\n//  ring tau = %s;", rString(currRing));
4526        rChangeCurrRing(oRing);
4527        G = idrMoveR(Gt, testring,currRing);
4528
4529        goto NEXT_VECTOR_FRACTAL;
4530#endif
4531
4532      FRACTAL_MSTDCC:
4533        //Print("\n//  apply BB-Alg in ring = %s;", rString(currRing));
4534        to=clock();
4535        G = MstdCC(Gt);
4536        xtextra=xtextra+clock()-to;
4537
4538        oRing = currRing;
4539
4540        // update the original target vector w.r.t. the current GB
4541        if(MivSame(Xivinput, Xivlp) == 1)
4542          if (rParameter(currRing) != NULL)
4543            DefRingParlp();
4544          else
4545            VMrDefaultlp();
4546        else
4547          if (rParameter(currRing) != NULL)
4548            DefRingPar(Xivinput);
4549          else
4550            VMrDefault(Xivinput);
4551
4552        testring = currRing;
4553        Gt = idrMoveR(G, oRing,currRing);
4554
4555        delete Xtau;
4556        Xtau = NewVectorlp(Gt);
4557
4558        rChangeCurrRing(oRing);
4559        G = idrMoveR(Gt, testring,currRing);
4560
4561        delete omega2;
4562        delete next_vect;
4563        delete altomega;
4564        /*
4565          Print("\n// Leaving the %d-th recursion with %d steps,", nlev,nwalks);
4566          Print(" ** Overflow_Error? (%d)", Overflow_Error);
4567        */
4568        if(Overflow_Error == TRUE)
4569          nnflow ++;
4570
4571        Overflow_Error = FALSE;
4572        return(G);
4573      }
4574    }
4575
4576    for(i=nV-1; i>=0; i--) {
4577      (*altomega)[i] = (*omega)[i];
4578      (*omega)[i] = (*next_vect)[i];
4579    }
4580    delete next_vect;
4581
4582    to=clock();
4583    /* Take the initial form of <G> w.r.t. omega */
4584    Gomega = MwalkInitialForm(G, omega);
4585    xtif=xtif+clock()-to;
4586
4587#ifndef  BUCHBERGER_ALG
4588    if(isNolVector(omega) == 0)
4589      hilb_func = hFirstSeries(Gomega,NULL,NULL,omega,currRing);
4590    else
4591      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
4592#endif // BUCHBERGER_ALG
4593
4594    if (rParameter(currRing) != NULL)
4595      DefRingPar(omega);
4596    else
4597      VMrDefault(omega);
4598
4599    Gomega1 = idrMoveR(Gomega, oRing,currRing);
4600
4601    /* Maximal recursion depth, to compute a red. GB */
4602    /* Fractal walk with the alternative recursion */
4603    /* alternative recursion */
4604    // if(nlev == nV || lengthpoly(Gomega1) == 0)
4605    if(nlev == Xnlev || lengthpoly(Gomega1) == 0)
4606      //if(nlev == nV) // blind recursion
4607    {
4608      /*
4609      if(Xnlev != nV)
4610      {
4611        Print("\n// ** Xnlev = %d", Xnlev);
4612        ivString(Xtau, "Xtau");
4613      }
4614      */
4615      to=clock();
4616#ifdef  BUCHBERGER_ALG
4617      Gresult = MstdhomCC(Gomega1);
4618#else
4619      Gresult =kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,omega);
4620      delete hilb_func;
4621#endif // BUCHBERGER_ALG
4622      xtstd=xtstd+clock()-to;
4623    }
4624    else {
4625      rChangeCurrRing(oRing);
4626      Gomega1 = idrMoveR(Gomega1, oRing,currRing);
4627      Gresult = rec_fractal_call(idCopy(Gomega1),nlev+1,omega);
4628    }
4629
4630    //convert a Groebner basis from a ring to another ring,
4631    new_ring = currRing;
4632
4633    rChangeCurrRing(oRing);
4634    Gresult1 = idrMoveR(Gresult, new_ring,currRing);
4635    Gomega2 = idrMoveR(Gomega1, new_ring,currRing);
4636
4637    to=clock();
4638    /* Lifting process */
4639    F = MLifttwoIdeal(Gomega2, Gresult1, G);
4640    xtlift=xtlift+clock()-to;
4641    idDelete(&Gresult1);
4642    idDelete(&Gomega2);
4643    idDelete(&G);
4644
4645    rChangeCurrRing(new_ring);
4646    F1 = idrMoveR(F, oRing,currRing);
4647
4648    to=clock();
4649    /* Interreduce G */
4650    G = kInterRedCC(F1, NULL);
4651    xtred=xtred+clock()-to;
4652    idDelete(&F1);
4653  }
4654}
4655
4656ideal Mfwalk(ideal G, intvec* ivstart, intvec* ivtarget)
4657{
4658  Set_Error(FALSE);
4659  Overflow_Error = FALSE;
4660  //Print("// pSetm_Error = (%d)", ErrorCheck());
4661  //Print("\n// ring ro = %s;", rString(currRing));
4662
4663  nnflow = 0;
4664  Xngleich = 0;
4665  Xcall = 0;
4666  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0;
4667  xftinput = clock();
4668
4669  ring  oldRing = currRing;
4670  int i, nV = currRing->N;
4671  XivNull = new intvec(nV);
4672  Xivinput = ivtarget;
4673  ngleich = 0;
4674  to=clock();
4675  ideal I = MstdCC(G);
4676  G = NULL;
4677  xftostd=clock()-to;
4678  Xsigma = ivstart;
4679
4680  Xnlev=nV;
4681
4682#ifdef FIRST_STEP_FRACTAL
4683  ideal Gw = MwalkInitialForm(I, ivstart);
4684  for(i=IDELEMS(Gw)-1; i>=0; i--)
4685  {
4686    if((Gw->m[i]!=NULL) /* len >=0 */
4687    && (Gw->m[i]->next!=NULL) /* len >=1 */
4688    && (Gw->m[i]->next->next!=NULL)) /* len >=2 */
4689    {
4690      intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
4691      intvec* Mdp;
4692
4693      if(MivSame(ivstart, iv_dp) != 1)
4694        Mdp = MivWeightOrderdp(ivstart);
4695      else
4696        Mdp = MivMatrixOrderdp(nV);
4697
4698      Xsigma = Mfpertvector(I, Mdp);
4699      Overflow_Error = FALSE;
4700
4701      delete Mdp;
4702      delete iv_dp;
4703      break;
4704    }
4705  }
4706  idDelete(&Gw);
4707#endif
4708
4709  ideal I1;
4710  intvec* Mlp;
4711  Xivlp = Mivlp(nV);
4712
4713  if(MivComp(ivtarget, Xivlp)  != 1)
4714  {
4715    if (rParameter(currRing) != NULL)
4716      DefRingPar(ivtarget);
4717    else
4718      VMrDefault(ivtarget);
4719
4720    I1 = idrMoveR(I, oldRing,currRing);
4721    Mlp = MivWeightOrderlp(ivtarget);
4722    Xtau = Mfpertvector(I1, Mlp);
4723  }
4724  else
4725  {
4726    if (rParameter(currRing) != NULL)
4727      DefRingParlp();
4728    else
4729      VMrDefaultlp();
4730
4731    I1 = idrMoveR(I, oldRing,currRing);
4732    Mlp =  MivMatrixOrderlp(nV);
4733    Xtau = Mfpertvector(I1, Mlp);
4734  }
4735  delete Mlp;
4736  Overflow_Error = FALSE;
4737
4738  //ivString(Xsigma, "Xsigma");
4739  //ivString(Xtau, "Xtau");
4740
4741  id_Delete(&I, oldRing);
4742  ring tRing = currRing;
4743
4744  if (rParameter(currRing) != NULL)
4745    DefRingPar(ivstart);
4746  else
4747    VMrDefault(ivstart);
4748
4749  I = idrMoveR(I1,tRing,currRing);
4750  to=clock();
4751  ideal J = MstdCC(I);
4752  idDelete(&I);
4753  xftostd=xftostd+clock()-to;
4754
4755  ideal resF;
4756  ring helpRing = currRing;
4757
4758  J = rec_fractal_call(J, 1, ivtarget);
4759
4760  rChangeCurrRing(oldRing);
4761  resF = idrMoveR(J, helpRing,currRing);
4762  idSkipZeroes(resF);
4763
4764  delete Xivlp;
4765  delete Xsigma;
4766  delete Xtau;
4767  delete XivNull;
4768
4769#ifdef TIME_TEST
4770  TimeStringFractal(xftinput, xftostd, xtif, xtstd, xtextra,
4771                    xtlift, xtred, xtnw);
4772
4773
4774  Print("\n// pSetm_Error = (%d)", ErrorCheck());
4775  Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
4776  Print("\n// the numbers of Overflow_Error (%d)", nnflow);
4777#endif
4778
4779  return(resF);
4780}
4781
4782/* Tran algorithm */
4783/* use kStd, if nP = 0, else call Ab_Rec_Pert (LastGB) */
4784ideal TranMImprovwalk(ideal G,intvec* curr_weight,intvec* target_tmp, int nP)
4785{
4786  clock_t mtim = clock();
4787  Set_Error(FALSE  );
4788  Overflow_Error =  FALSE;
4789  //Print("// pSetm_Error = (%d)", ErrorCheck());
4790  //Print("\n// ring ro = %s;", rString(currRing));
4791
4792  clock_t tostd, tif=0, tstd=0, tlift=0, tred=0, tnw=0, textra=0;
4793  clock_t tinput = clock();
4794
4795  int nsteppert=0, i, nV = currRing->N, nwalk=0, npert_tmp=0;
4796  int *npert=(int*)omAlloc(2*nV*sizeof(int));
4797  ideal Gomega, M,F,  G1, Gomega1, Gomega2, M1, F1;
4798  ring endRing, newRing, oldRing, lpRing;
4799  intvec* next_weight;
4800  intvec* ivNull = new intvec(nV); //define (0,...,0)
4801  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
4802  intvec* iv_lp = Mivlp(nV); //define (1,0,...,0)
4803  ideal H0, H1, H2, Glp;
4804  int nGB, endwalks = 0,  nwalkpert=0,  npertstep=0;
4805  intvec* Mlp =  MivMatrixOrderlp(nV);
4806  intvec* vector_tmp = new intvec(nV);
4807  intvec* hilb_func;
4808
4809  /* to avoid (1,0,...,0) as the target vector */
4810  intvec* last_omega = new intvec(nV);
4811  for(i=nV-1; i>0; i--)
4812    (*last_omega)[i] = 1;
4813  (*last_omega)[0] = 10000;
4814
4815  //  intvec* extra_curr_weight = new intvec(nV);
4816  intvec* target_weight = new intvec(nV);
4817  for(i=nV-1; i>=0; i--)
4818    (*target_weight)[i] = (*target_tmp)[i];
4819
4820  ring XXRing = currRing;
4821  newRing = currRing;
4822
4823  to=clock();
4824  /* compute a red. GB w.r.t. the help ring */
4825  if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) = "dp"
4826    G = MstdCC(G);
4827  else
4828  {
4829    //rOrdStr(currRing) = (a(.c_w..),lp,C)
4830    if (rParameter(currRing) != NULL)
4831      DefRingPar(curr_weight);
4832    else
4833      VMrDefault(curr_weight);
4834    G = idrMoveR(G, XXRing,currRing);
4835    G = MstdCC(G);
4836  }
4837  tostd=clock()-to;
4838
4839#ifdef REPRESENTATION_OF_SIGMA
4840  ideal Gw = MwalkInitialForm(G, curr_weight);
4841
4842  if(islengthpoly2(Gw)==1)
4843  {
4844    intvec* MDp;
4845    if(MivComp(curr_weight, iv_dp) == 1)
4846      MDp = MatrixOrderdp(nV); //MivWeightOrderlp(iv_dp);
4847    else
4848      MDp = MivWeightOrderlp(curr_weight);
4849
4850    curr_weight = RepresentationMatrix_Dp(G, MDp);
4851
4852    delete MDp;
4853
4854    ring exring = currRing;
4855
4856    if (rParameter(currRing) != NULL)
4857      DefRingPar(curr_weight);
4858    else
4859      VMrDefault(curr_weight);
4860    to=clock();
4861    Gw = idrMoveR(G, exring,currRing);
4862    G = MstdCC(Gw);
4863    Gw = NULL;
4864    tostd=tostd+clock()-to;
4865    //ivString(curr_weight,"rep. sigma");
4866    goto COMPUTE_NEW_VECTOR;
4867  }
4868
4869  idDelete(&Gw);
4870  delete iv_dp;
4871#endif
4872
4873
4874  while(1)
4875  {
4876    to=clock();
4877    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
4878    Gomega = MwalkInitialForm(G, curr_weight);
4879    tif=tif+clock()-to;
4880
4881#ifndef  BUCHBERGER_ALG
4882    if(isNolVector(curr_weight) == 0)
4883      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
4884    else
4885      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
4886#endif // BUCHBERGER_ALG
4887
4888    oldRing = currRing;
4889
4890    /* define a new ring that its ordering is "(a(curr_weight),lp) */
4891    if (rParameter(currRing) != NULL)
4892      DefRingPar(curr_weight);
4893    else
4894      VMrDefault(curr_weight);
4895
4896    newRing = currRing;
4897    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
4898
4899    to=clock();
4900    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
4901#ifdef  BUCHBERGER_ALG
4902    M = MstdhomCC(Gomega1);
4903#else
4904    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
4905    delete hilb_func;
4906#endif // BUCHBERGER_ALG
4907    tstd=tstd+clock()-to;
4908
4909    /* change the ring to oldRing */
4910    rChangeCurrRing(oldRing);
4911    M1 =  idrMoveR(M, newRing,currRing);
4912    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
4913
4914    to=clock();
4915    /* compute a representation of the generators of submod (M)
4916       with respect to those of mod (Gomega).
4917       Gomega is a reduced Groebner basis w.r.t. the current ring */
4918    F = MLifttwoIdeal(Gomega2, M1, G);
4919    tlift=tlift+clock()-to;
4920
4921    idDelete(&M1);
4922    idDelete(&Gomega2);
4923    idDelete(&G);
4924
4925    /* change the ring to newRing */
4926    rChangeCurrRing(newRing);
4927    F1 = idrMoveR(F, oldRing,currRing);
4928
4929    to=clock();
4930    /* reduce the Groebner basis <G> w.r.t. new ring */
4931    G = kInterRedCC(F1, NULL);
4932    tred=tred+clock()-to;
4933    idDelete(&F1);
4934
4935
4936  COMPUTE_NEW_VECTOR:
4937    newRing = currRing;
4938    nwalk++;
4939    nwalkpert++;
4940    to=clock();
4941    /* compute a next weight vector */
4942    next_weight = MwalkNextWeightCC(curr_weight,target_weight, G);
4943    tnw=tnw+clock()-to;
4944#ifdef PRINT_VECTORS
4945    MivString(curr_weight, target_weight, next_weight);
4946#endif
4947
4948
4949    /* check whether the computed intermediate weight vector in
4950       the correct cone is, since sometimes it is very big e.g. s7, cyc7.
4951       If it is NOT in the cone, then one has directly to compute
4952       a reduced Groebner basis with respect to the lexicographic ordering
4953       for the known Groebner basis that it is computed in the last step.
4954    */
4955    //if(test_w_in_ConeCC(G, next_weight) != 1)
4956    if(Overflow_Error == TRUE)
4957    {
4958    OMEGA_OVERFLOW_TRAN_NEW:
4959      //Print("\n//  takes %d steps!", nwalk-1);
4960      //Print("\n//ring lastRing = %s;", rString(currRing));
4961#ifdef TEST_OVERFLOW
4962      goto  BE_FINISH;
4963#endif
4964
4965#ifdef CHECK_IDEAL_MWALK
4966      idElements(G, "G");
4967      //headidString(G, "G");
4968#endif
4969
4970      if(MivSame(target_tmp, iv_lp) == 1)
4971        if (rParameter(currRing) != NULL)
4972          DefRingParlp();
4973        else
4974          VMrDefaultlp();
4975      else
4976        if (rParameter(currRing) != NULL)
4977          DefRingPar(target_tmp);
4978        else
4979          VMrDefault(target_tmp);
4980
4981      lpRing = currRing;
4982      G1 = idrMoveR(G, newRing,currRing);
4983
4984      to=clock();
4985      /*apply kStd or LastGB to compute  a lex. red. Groebner basis of <G>*/
4986      if(nP == 0 || MivSame(target_tmp, iv_lp) == 0){
4987        //Print("\n\n// calls \"std in ring r_%d = %s;", nwalk, rString(currRing));
4988        G = MstdCC(G1);//no result for qnt1
4989      }
4990      else {
4991        rChangeCurrRing(newRing);
4992        G1 = idrMoveR(G1, lpRing,currRing);
4993
4994        //Print("\n\n// calls \"LastGB\" (%d) to compute a GB", nV-1);
4995        G = LastGB(G1, curr_weight, nV-1); //no result for kats7
4996
4997        rChangeCurrRing(lpRing);
4998        G = idrMoveR(G, newRing,currRing);
4999      }
5000      textra=clock()-to;
5001      npert[endwalks]=nwalk-npert_tmp;
5002      npert_tmp = nwalk;
5003      endwalks ++;
5004      break;
5005    }
5006
5007    /* check whether the computed Groebner basis a really Groebner basis is.
5008       if no, we perturb the target vector with the maximal "perturbation"
5009       degree.*/
5010    if(MivComp(next_weight, target_weight) == 1 ||
5011       MivComp(next_weight, curr_weight) == 1 )
5012    {
5013      //Print("\n//ring r_%d = %s;", nwalk, rString(currRing));
5014
5015
5016      //compute the number of perturbations and its step
5017      npert[endwalks]=nwalk-npert_tmp;
5018      npert_tmp = nwalk;
5019
5020      endwalks ++;
5021
5022      /*it is very important if the walk only uses one step, e.g. Fate, liu*/
5023      if(endwalks == 1 && MivComp(next_weight, curr_weight) == 1){
5024        rChangeCurrRing(XXRing);
5025        G = idrMoveR(G, newRing,currRing);
5026        goto FINISH;
5027      }
5028      H0 = id_Head(G,currRing);
5029
5030      if(MivSame(target_tmp, iv_lp) == 1)
5031        if (rParameter(currRing) != NULL)
5032          DefRingParlp();
5033        else
5034          VMrDefaultlp();
5035      else
5036        if (rParameter(currRing) != NULL)
5037          DefRingPar(target_tmp);
5038        else
5039          VMrDefault(target_tmp);
5040
5041      lpRing = currRing;
5042      Glp = idrMoveR(G, newRing,currRing);
5043      H2 = idrMoveR(H0, newRing,currRing);
5044
5045      /* Apply Lemma 2.2 in Collart et. al (1997) to check whether
5046         cone(k-1) is equal to cone(k) */
5047      nGB = 1;
5048      for(i=IDELEMS(Glp)-1; i>=0; i--)
5049      {
5050        poly t;
5051        if((t=pSub(pHead(Glp->m[i]), pCopy(H2->m[i]))) != NULL)
5052        {
5053          pDelete(&t);
5054          idDelete(&H2);//5.5.02
5055          nGB = 0; //i.e. Glp is no reduced Groebner basis
5056          break;
5057        }
5058        pDelete(&t);
5059      }
5060
5061      idDelete(&H2);//5.5.02
5062
5063      if(nGB == 1)
5064      {
5065        G = Glp;
5066        Glp = NULL;
5067        break;
5068      }
5069
5070       /* perturb the target weight vector, if the vector target_tmp
5071          stays in many cones */
5072      poly p;
5073      BOOLEAN plength3 = FALSE;
5074      for(i=IDELEMS(Glp)-1; i>=0; i--)
5075      {
5076        p = MpolyInitialForm(Glp->m[i], target_tmp);
5077        if(p->next != NULL &&
5078           p->next->next != NULL &&
5079           p->next->next->next != NULL)
5080        {
5081          Overflow_Error = FALSE;
5082
5083          for(i=0; i<nV; i++)
5084            (*vector_tmp)[i] = (*target_weight)[i];
5085
5086          delete target_weight;
5087          target_weight = MPertVectors(Glp, Mlp, nV);
5088
5089          if(MivComp(vector_tmp, target_weight)==1)
5090          {
5091            //PrintS("\n// The old and new representaion vector are the same!!");
5092            G = Glp;
5093            newRing = currRing;
5094            goto OMEGA_OVERFLOW_TRAN_NEW;
5095           }
5096
5097          if(Overflow_Error == TRUE)
5098          {
5099            rChangeCurrRing(newRing);
5100            G = idrMoveR(Glp, lpRing,currRing);
5101            goto OMEGA_OVERFLOW_TRAN_NEW;
5102          }
5103
5104          plength3 = TRUE;
5105          pDelete(&p);
5106          break;
5107        }
5108        pDelete(&p);
5109      }
5110
5111      if(plength3 == FALSE)
5112      {
5113        rChangeCurrRing(newRing);
5114        G = idrMoveR(Glp, lpRing,currRing);
5115        goto TRAN_LIFTING;
5116      }
5117
5118
5119      npertstep = nwalk;
5120      nwalkpert = 1;
5121      nsteppert ++;
5122
5123      /*
5124      Print("\n// Subroutine needs (%d) steps.", nwalk);
5125      idElements(Glp, "last G in walk:");
5126      PrintS("\n// ****************************************");
5127      Print("\n// Perturb the original target vector (%d): ", nsteppert);
5128      ivString(target_weight, "new target");
5129      PrintS("\n// ****************************************\n");
5130      */
5131      rChangeCurrRing(newRing);
5132      G = idrMoveR(Glp, lpRing,currRing);
5133
5134      delete next_weight;
5135
5136      //Print("\n// ring rNEW = %s;", rString(currRing));
5137      goto COMPUTE_NEW_VECTOR;
5138    }
5139
5140  TRAN_LIFTING:
5141    for(i=nV-1; i>=0; i--)
5142      (*curr_weight)[i] = (*next_weight)[i];
5143
5144    delete next_weight;
5145  }//while
5146
5147 BE_FINISH:
5148  rChangeCurrRing(XXRing);
5149  G = idrMoveR(G, lpRing,currRing);
5150
5151 FINISH:
5152  delete ivNull;
5153  delete next_weight;
5154  delete iv_lp;
5155  omFree(npert);
5156
5157#ifdef TIME_TEST
5158  Print("\n// Computation took %d steps and %.2f sec",
5159        nwalk, ((double) (clock()-mtim)/1000000));
5160
5161  TimeStringFractal(tinput, tostd, tif, tstd, textra, tlift, tred, tnw);
5162
5163  Print("\n// pSetm_Error = (%d)", ErrorCheck());
5164  Print("\n// Overflow_Error? (%d)\n", Overflow_Error);
5165#endif
5166
5167  return(G);
5168}
5169
5170
5171/* compute the reduced Grï¿œbner basis of an ideal <Go> w.r.t. lp */
5172static ideal Mpwalk_MAltwalk1(ideal Go, intvec* curr_weight, int tp_deg)
5173{
5174  Overflow_Error = FALSE;
5175  BOOLEAN nOverflow_Error = FALSE;
5176  clock_t tproc=0;
5177  clock_t tinput=clock();
5178  int i, nV = currRing->N;
5179  int nwalk=0, endwalks=0, ntestwinC=1;
5180  int tp_deg_tmp = tp_deg;
5181  ideal Gomega, M, F, G, M1, F1, Gomega1, Gomega2, G1;
5182  ring endRing, newRing, oldRing, TargetRing;
5183  intvec* next_weight;
5184  intvec* ivNull = new intvec(nV);
5185  intvec* extra_curr_weight = new intvec(nV);
5186
5187  ring YXXRing = currRing;
5188
5189  intvec* iv_M_dpp = MivMatrixOrderlp(nV);
5190  intvec* target_weight;// = Mivlp(nV);
5191  ideal ssG;
5192
5193  /* perturb the target vector */
5194  while(1)
5195  {
5196    if(Overflow_Error == FALSE)
5197    {
5198      if (rParameter(currRing) != NULL)
5199        DefRingParlp();
5200      else
5201        VMrDefaultlp();
5202
5203      TargetRing = currRing;
5204      ssG = idrMoveR(Go,YXXRing,currRing);
5205    }
5206    Overflow_Error = FALSE;
5207    if(tp_deg != 1)
5208      target_weight = MPertVectors(ssG, iv_M_dpp, tp_deg);
5209    else
5210    {
5211      target_weight = Mivlp(nV);
5212      break;
5213    }
5214    if(Overflow_Error == FALSE)
5215      break;
5216
5217    Overflow_Error = TRUE;
5218    tp_deg --;
5219  }
5220  if(tp_deg != tp_deg_tmp)
5221  {
5222    Overflow_Error = TRUE;
5223    nOverflow_Error = TRUE;
5224  }
5225
5226  //  Print("\n// tp_deg = %d", tp_deg);
5227  // ivString(target_weight, "pert target");
5228
5229  delete iv_M_dpp;
5230  intvec* hilb_func;
5231
5232  /* to avoid (1,0,...,0) as the target vector */
5233  intvec* last_omega = new intvec(nV);
5234  for(i=nV-1; i>0; i--)
5235    (*last_omega)[i] = 1;
5236  (*last_omega)[0] = 10000;
5237
5238  rChangeCurrRing(YXXRing);
5239  G = idrMoveR(ssG, TargetRing,currRing);
5240
5241  while(1)
5242  {
5243    nwalk ++;
5244    nstep ++;
5245
5246    if(nwalk==1)
5247      goto FIRST_STEP;
5248
5249    to=clock();
5250    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
5251    Gomega = MwalkInitialForm(G, curr_weight);
5252    xtif=xtif+clock()-to;
5253#if 0
5254    if(Overflow_Error == TRUE)
5255    {
5256      for(i=nV-1; i>=0; i--)
5257        (*curr_weight)[i] = (*extra_curr_weight)[i];
5258      delete extra_curr_weight;
5259      goto LASTGB_ALT1;
5260    }
5261#endif
5262#ifndef  BUCHBERGER_ALG
5263    if(isNolVector(curr_weight) == 0)
5264      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
5265    else
5266      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
5267#endif // BUCHBERGER_ALG
5268
5269    oldRing = currRing;
5270
5271    /* define a new ring that its ordering is "(a(curr_weight),lp) */
5272    if (rParameter(currRing) != NULL)
5273      DefRingPar(curr_weight);
5274    else
5275      VMrDefault(curr_weight);
5276
5277    newRing = currRing;
5278    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
5279
5280#ifdef ENDWALKS
5281    if(endwalks == 1)
5282    {
5283      Print("\n//  it is  %d-th step!!", nwalk);
5284      idElements(Gomega1, "Gw");
5285      PrintS("\n//  compute a rGB of Gw:");
5286    }
5287#endif
5288
5289    to=clock();
5290    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
5291#ifdef  BUCHBERGER_ALG
5292    M = MstdhomCC(Gomega1);
5293#else
5294    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
5295    delete hilb_func;
5296#endif // BUCHBERGER_ALG
5297    xtstd=xtstd+clock()-to;
5298
5299    /* change the ring to oldRing */
5300    rChangeCurrRing(oldRing);
5301    M1 =  idrMoveR(M, newRing,currRing);
5302    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
5303    to=clock();
5304
5305    // if(endwalks == 1) PrintS("\n//  Lifting is still working:");
5306
5307    /* compute a reduced Groebner basis of <G> w.r.t. "newRing" by the
5308     lifting process */
5309    F = MLifttwoIdeal(Gomega2, M1, G);
5310    xtlift=xtlift+clock()-to;
5311
5312    idDelete(&M1);
5313    idDelete(&Gomega2);
5314    idDelete(&G);
5315
5316    /* change the ring to newRing */
5317    rChangeCurrRing(newRing);
5318    F1 = idrMoveR(F, oldRing,currRing);
5319    to=clock();
5320    //if(endwalks == 1) PrintS("\n//  InterRed is still working:");
5321    /* reduce the Groebner basis <G> w.r.t. the new ring */
5322    G = kInterRedCC(F1, NULL);
5323    xtred=xtred+clock()-to;
5324    idDelete(&F1);
5325
5326    if(endwalks == 1)
5327      break;
5328
5329  FIRST_STEP:
5330    Overflow_Error=FALSE;
5331    to=clock();
5332    /* compute a next weight vector */
5333    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
5334    xtnw=xtnw+clock()-to;
5335#ifdef PRINT_VECTORS
5336    MivString(curr_weight, target_weight, next_weight);
5337#endif
5338
5339    if(Overflow_Error == TRUE)
5340    {
5341      delete next_weight;
5342
5343      LASTGB_ALT1:
5344      if(tp_deg > 1){
5345        nOverflow_Error = Overflow_Error;
5346        tproc = tproc+clock()-tinput;
5347        /*
5348          Print("\n// A subroutine takes %d steps and calls \"Mpwalk\" (1,%d):",
5349          nwalk, tp_deg-1);
5350        */
5351        G1 = Mpwalk_MAltwalk1(G, curr_weight, tp_deg-1);
5352        goto MPW_Finish;
5353      }
5354      else {
5355        newRing = currRing;
5356        ntestwinC = 0;
5357        break;
5358      }
5359    }
5360
5361    if(MivComp(next_weight, ivNull) == 1)
5362    {
5363      newRing = currRing;
5364      delete next_weight;
5365      break;
5366    }
5367    if(MivComp(next_weight, target_weight) == 1)
5368      endwalks = 1;
5369
5370    for(i=nV-1; i>=0; i--)
5371    {
5372      //(*extra_curr_weight)[i] = (*curr_weight)[i];
5373      (*curr_weight)[i] = (*next_weight)[i];
5374    }
5375    delete next_weight;
5376  }//while
5377
5378  /* check wheather the pertubed target vector is correct */
5379
5380  //define and execute ring with lex. order
5381  if (rParameter(currRing) != NULL)
5382    DefRingParlp();
5383  else
5384    VMrDefaultlp();
5385
5386  G1 = idrMoveR(G, newRing,currRing);
5387
5388  if( test_w_in_ConeCC(G1, target_weight) != 1 || ntestwinC == 0)
5389  {
5390    PrintS("\n// The perturbed target vector doesn't STAY in the correct cone!!");
5391    if(tp_deg == 1){
5392      //Print("\n// subroutine takes %d steps and applys \"std\"", nwalk);
5393      to=clock();
5394      ideal G2 = MstdCC(G1);
5395      xtextra=xtextra+clock()-to;
5396      idDelete(&G1);
5397      G1 = G2;
5398      G2 = NULL;
5399    }
5400    else {
5401      nOverflow_Error = Overflow_Error;
5402      tproc = tproc+clock()-tinput;
5403      /*
5404        Print("\n// B subroutine takes %d steps and calls \"Mpwalk\" (1,%d) :",
5405            nwalk,  tp_deg-1);
5406      */
5407      G1 = Mpwalk_MAltwalk1(G1, curr_weight, tp_deg-1);
5408    }
5409  }
5410
5411 MPW_Finish:
5412  newRing = currRing;
5413  rChangeCurrRing(YXXRing);
5414  ideal result = idrMoveR(G1, newRing,currRing);
5415
5416  delete ivNull;
5417  delete target_weight;
5418
5419  /*
5420  Print("\n// \"Mpwalk\" (1,%d) took %d steps and %.2f sec. Overflow_Error (%d)", tp_deg,
5421        nwalk, ((double) clock()-tinput)/1000000, nOverflow_Error);
5422  */
5423
5424  return(result);
5425}
5426
5427/* August 2003 */
5428ideal MAltwalk1(ideal Go, int op_deg, int tp_deg, intvec* curr_weight,
5429                intvec* target_weight)
5430{
5431  Set_Error(FALSE  );
5432  Overflow_Error = FALSE;
5433  BOOLEAN nOverflow_Error = FALSE;
5434  // Print("// pSetm_Error = (%d)", ErrorCheck());
5435
5436  xtif=0; xtstd=0; xtlift=0; xtred=0; xtnw=0; xtextra=0;
5437  xftinput = clock();
5438  clock_t tostd, tproc;
5439
5440  nstep = 0;
5441  int i, nV = currRing->N;
5442  int nwalk=0, endwalks=0;
5443  int op_tmp = op_deg;
5444  ideal Gomega, M, F, G, G1, Gomega1, Gomega2, M1, F1;
5445  ring endRing, newRing, oldRing, TargetRing;
5446  intvec* next_weight;
5447  intvec* iv_M_dp;
5448  intvec* ivNull = new intvec(nV);
5449  intvec* iv_dp = MivUnit(nV);// define (1,1,...,1)
5450  intvec* exivlp = Mivlp(nV);
5451  intvec* extra_curr_weight = new intvec(nV);
5452  intvec* hilb_func;
5453
5454  intvec* cw_tmp = curr_weight;
5455
5456  /* to avoid (1,0,...,0) as the target vector */
5457  intvec* last_omega = new intvec(nV);
5458  for(i=nV-1; i>0; i--)
5459    (*last_omega)[i] = 1;
5460  (*last_omega)[0] = 10000;
5461
5462  ring XXRing = currRing;
5463
5464  to=clock();
5465  /* compute a pertubed weight vector of the original weight vector.
5466     The perturbation degree is recursive decrease until that vector
5467     stays inn the correct cone. */
5468  while(1)
5469  {
5470    if(Overflow_Error == FALSE)
5471    {
5472      if(MivComp(curr_weight, iv_dp) == 1) //rOrdStr(currRing) = "dp"
5473        if(op_tmp == op_deg) {
5474          G = MstdCC(Go);
5475          if(op_deg != 1)
5476            iv_M_dp = MivMatrixOrderdp(nV);
5477        }
5478    }
5479    else
5480    {
5481      if(op_tmp == op_deg) {
5482        //rOrdStr(currRing) = (a(...),lp,C)
5483        if (rParameter(currRing) != NULL)
5484          DefRingPar(cw_tmp);
5485        else
5486          VMrDefault(cw_tmp);
5487
5488        G = idrMoveR(Go, XXRing,currRing);
5489        G = MstdCC(G);
5490        if(op_deg != 1)
5491          iv_M_dp = MivMatrixOrder(cw_tmp);
5492      }
5493    }
5494    Overflow_Error = FALSE;
5495    if(op_deg != 1)
5496      curr_weight = MPertVectors(G, iv_M_dp, op_deg);
5497    else {
5498      curr_weight =  cw_tmp;
5499      break;
5500    }
5501    if(Overflow_Error == FALSE)
5502      break;
5503
5504    Overflow_Error = TRUE;
5505    op_deg --;
5506  }
5507  tostd=clock()-to;
5508
5509  if(op_tmp != 1 )
5510    delete iv_M_dp;
5511  delete iv_dp;
5512
5513  if(currRing->order[0] == ringorder_a)
5514    goto NEXT_VECTOR;
5515
5516  while(1)
5517  {
5518    nwalk ++;
5519    nstep ++;
5520
5521    to = clock();
5522    /* compute an initial form ideal of <G> w.r.t. "curr_vector" */
5523    Gomega = MwalkInitialForm(G, curr_weight);
5524    xtif=xtif+clock()-to;
5525#if 0
5526    if(Overflow_Error == TRUE)
5527    {
5528      for(i=nV-1; i>=0; i--)
5529        (*curr_weight)[i] = (*extra_curr_weight)[i];
5530      delete extra_curr_weight;
5531
5532      newRing = currRing;
5533      goto MSTD_ALT1;
5534    }
5535#endif
5536#ifndef  BUCHBERGER_ALG
5537    if(isNolVector(curr_weight) == 0)
5538      hilb_func = hFirstSeries(Gomega,NULL,NULL,curr_weight,currRing);
5539    else
5540      hilb_func = hFirstSeries(Gomega,NULL,NULL,last_omega,currRing);
5541#endif // BUCHBERGER_ALG
5542
5543    oldRing = currRing;
5544
5545    /* define a new ring which ordering is "(a(curr_weight),lp) */
5546    if (rParameter(currRing) != NULL)
5547      DefRingPar(curr_weight);
5548    else
5549      VMrDefault(curr_weight);
5550
5551    newRing = currRing;
5552    Gomega1 = idrMoveR(Gomega, oldRing,currRing);
5553
5554    to=clock();
5555    /* compute a reduced Groebner basis of <Gomega> w.r.t. "newRing" */
5556#ifdef  BUCHBERGER_ALG
5557    M = MstdhomCC(Gomega1);
5558#else
5559    M=kStd(Gomega1,NULL,isHomog,NULL,hilb_func,0,NULL,curr_weight);
5560    delete hilb_func;
5561#endif // BUCHBERGER_ALG
5562    xtstd=xtstd+clock()-to;
5563
5564    /* change the ring to oldRing */
5565    rChangeCurrRing(oldRing);
5566    M1 =  idrMoveR(M, newRing,currRing);
5567    Gomega2 =  idrMoveR(Gomega1, newRing,currRing);
5568
5569    to=clock();
5570    /* compute a reduced Groebner basis of <G> w.r.t. "newRing" by the
5571     lifting process */
5572    F = MLifttwoIdeal(Gomega2, M1, G);
5573    xtlift=xtlift+clock()-to;
5574
5575    idDelete(&M1);
5576    idDelete(&Gomega2);
5577    idDelete(&G);
5578
5579   /* change the ring to newRing */
5580    rChangeCurrRing(newRing);
5581    F1 = idrMoveR(F, oldRing,currRing);
5582
5583    to=clock();
5584    /* reduce the Groebner basis <G> w.r.t. new ring */
5585    G = kInterRedCC(F1, NULL);
5586    xtred=xtred+clock()-to;
5587    idDelete(&F1);
5588
5589    if(endwalks == 1)
5590      break;
5591
5592  NEXT_VECTOR:
5593    to=clock();
5594    /* compute a next weight vector */
5595    next_weight = MkInterRedNextWeight(curr_weight,target_weight, G);
5596    xtnw=xtnw+clock()-to;
5597#ifdef PRINT_VECTORS
5598    MivString(curr_weight, target_weight, next_weight);
5599#endif
5600
5601    if(Overflow_Error == TRUE)
5602    {
5603      newRing = currRing;
5604
5605      if (rParameter(currRing) != NULL)
5606        DefRingPar(target_weight);
5607      else
5608        VMrDefault(target_weight);
5609
5610      F1 = idrMoveR(G, newRing,currRing);
5611      G = MstdCC(F1);
5612      idDelete(&F1);
5613      newRing = currRing;
5614      break; //for while
5615    }
5616
5617
5618    /* G is the wanted Groebner basis if next_weight == curr_weight */
5619    if(MivComp(next_weight, ivNull) == 1)
5620    {
5621      newRing = currRing;
5622      delete next_weight;
5623      break; //for while
5624    }
5625
5626    if(MivComp(next_weight, target_weight) == 1)
5627    {
5628      if(tp_deg == 1 || MivSame(target_weight, exivlp) == 0)
5629        endwalks = 1;
5630      else
5631      {
5632        MSTD_ALT1:
5633        nOverflow_Error = Overflow_Error;
5634        tproc = clock()-xftinput;
5635        /*
5636        Print("\n//  main routine takes %d steps and calls \"Mpwalk\" (1,%d):",
5637              nwalk,  tp_deg);
5638        */
5639        // compute the red. GB of <G> w.r.t. the lex order by
5640        // the "recursive-modified" perturbation walk alg (1,tp_deg)
5641        G = Mpwalk_MAltwalk1(G, curr_weight, tp_deg);
5642        delete next_weight;
5643        break; // for while
5644      }
5645    }
5646
5647    /* 06.11.01 NOT Changed, to free memory*/
5648    for(i=nV-1; i>=0; i--)
5649    {
5650      //(*extra_curr_weight)[i] = (*curr_weight)[i];
5651      (*curr_weight)[i] = (*next_weight)[i];
5652    }
5653    delete next_weight;
5654  }//while
5655
5656  rChangeCurrRing(XXRing);
5657  ideal result = idrMoveR(G, newRing,currRing);
5658  id_Delete(&G, newRing);
5659
5660  delete ivNull;
5661  if(op_deg != 1 )
5662    delete curr_weight;
5663
5664  delete exivlp;
5665#ifdef TIME_TEST
5666
5667  Print("\n// \"Main procedure\"  took %d steps, %.2f sec. and Overflow_Error(%d)",
5668        nwalk, ((double) tproc)/1000000, nOverflow_Error);
5669
5670  TimeStringFractal(xftinput, tostd, xtif, xtstd,xtextra, xtlift, xtred, xtnw);
5671
5672  Print("\n// pSetm_Error = (%d)", ErrorCheck());
5673  Print("\n// Overflow_Error? (%d)", Overflow_Error);
5674  Print("\n// Awalk1 took %d steps.\n", nstep);
5675#endif
5676  return(result);
5677}
5678
Note: See TracBrowser for help on using the repository browser.