source: git/Singular/walk.cc @ 6ce030f

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