source: git/Singular/walk.cc @ 667ba1

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