source: git/Singular/walk.cc @ f23590

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