source: git/Singular/walk.cc @ 993ae2

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