source: git/Singular/walk.cc @ ed8227

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