source: git/kernel/walkMain.cc @ 88479ff

spielwiese
Last change on this file since 88479ff was 57592b3, checked in by Hans Schoenemann <hannes@…>, 13 years ago
fix includes for int64vec.h
  • Property mode set to 100644
File size: 17.5 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5/*
6 * ABSTRACT: fractal walk stuff
7 *
8*/
9#include <string.h>
10#include <kernel/mod2.h>
11#include <misc/options.h>
12#include <misc/intvec.h>
13#include <misc/int64vec.h>
14#include <polys/polys.h>
15#include <kernel/ideals.h>
16#include <polys/monomials/ring.h>
17#include <kernel/walkMain.h>
18#include <kernel/walkSupport.h>
19#include <polys/prCopy.h>
20#include <kernel/kstd1.h>
21#include <polys/matpol.h>
22#include <polys/monomials/ring.h>
23
24
25///////////////////////////////////////////////////////////////////
26//Groebner Walk and Fractal Walk
27///////////////////////////////////////////////////////////////////
28//v1.3 2004-11-15
29///////////////////////////////////////////////////////////////////
30//implemented by Henrik Strohmayer
31///////////////////////////////////////////////////////////////////
32//in C++:
33///////////////////////////////////////////////////////////////////
34
35
36int overflow_error; //global variable
37
38/*
39overflow_error table
40 1: Miv64DotProduct mult
41 2: Miv64DotProduct add
42 3: gett64 zaehler mult
43 4: gett64 zaehler add
44 5: gett64 nenner mult
45 6: gett64 nenner add
46 7: nextw64 mult a
47 8: nextw64 mult b
48 9: nextw64 add
4910: getinveps64 mult
5011: getinveps64 add
5112: getTaun64 mult
5213: getTaun64 add
53*/
54
55///////////////////////////////////////////////////////////////////
56//fistWalkStep64
57///////////////////////////////////////////////////////////////////
58//Description: adapts currRing for the algorithm and moves G into
59//it, calculating the appropriate GB if currw not is inside of
60//the current groebner cone
61///////////////////////////////////////////////////////////////////
62//Assumes: that the option redSB is turned off
63///////////////////////////////////////////////////////////////////
64//Uses: init64, rCopy0AndAddA, rComplete, rChangeCurrRing, matIdLift,
65//mpMult, idInterRed, idStd, idCopy, idrMoveR,currwOnBorder64
66///////////////////////////////////////////////////////////////////
67
68WalkState firstWalkStep64(ideal & G,int64vec* currw64, ring destRing){
69  WalkState state=WalkOk;
70  /* OLDRING **************************************************** */
71  ideal nextG;
72  BITSET optionState;
73
74  if (currwOnBorder64(G,currw64))
75  {
76    ideal Gw=init64(G,currw64);
77    ring oldRing=currRing;
78  /* NEWRING **************************************************** */
79    ring rnew=rCopy0AndAddA(destRing,currw64);
80    rComplete(rnew);
81    rChangeCurrRing(rnew);
82
83    ideal newGw=idrMoveR(Gw, oldRing, rnew);
84
85
86 //HIER GEANDERT
87  matrix L=mpNew(1,1);
88  idLiftStd(newGw,&L);
89
90//       //turn off bucket representation of polynomials and on redSB
91//     optionState=test;
92//     //test|=Sy_bit(OPT_NOT_BUCKETS);
93//     test|=Sy_bit(OPT_REDSB);
94
95//     ideal newStdGw=idStd(newGw);
96
97//    //turn on bucket representation of polynomials and off redSB
98//    test=optionState;
99
100//     matrix L=matIdLift(newGw,newStdGw);
101//     idDelete(&newStdGw);
102
103    idDelete(&newGw);
104
105    nextG=idrMoveR(G,oldRing, rnew); idTest(nextG);
106
107    matrix nextGmat=(matrix)nextG;
108
109    matrix resMat=mp_Mult(nextGmat,L,rnew);
110    idDelete((ideal *)&nextGmat);
111    idDelete((ideal *)&L);
112
113    nextG=(ideal)resMat;
114
115    optionState=test;
116    test|=Sy_bit(OPT_REDSB);
117    nextG = idInterRed(nextG);
118    test=optionState;
119  }
120  else
121  {
122    ring oldRing=currRing;
123    ring rnew=rCopy0AndAddA(destRing,currw64);
124    rComplete(rnew);
125    rChangeCurrRing(rnew);
126    nextG=idrMoveR(G,oldRing, rnew);
127  }
128
129  G=nextG;
130  return(state);
131}
132
133///////////////////////////////////////////////////////////////////
134
135
136///////////////////////////////////////////////////////////////////
137//walkStep64
138///////////////////////////////////////////////////////////////////
139//Description: one step in the groebner walk
140///////////////////////////////////////////////////////////////////
141//Assumes: that the option redSB is turned off
142///////////////////////////////////////////////////////////////////
143//Uses: init, rCopyAndChangeA, matIdLift, mpMult,
144//idInterRed, mStd, idCopy, idrMoveR
145///////////////////////////////////////////////////////////////////
146
147WalkState walkStep64(ideal & G,int64vec* currw64, int step){
148  WalkState state=WalkOk;
149  BITSET optionState;
150
151/* OLDRING ****************************************************** */
152  ideal Gw=init64(G,currw64);
153
154  ring oldRing=currRing;
155
156/* NEWRING ****************************************************** */
157  rCopyAndChangeA(currw64);
158
159  ideal newGw=idrMoveR(Gw, oldRing,currRing);
160
161  //HIER GEANDERT
162  matrix L=mpNew(1,1);
163  ideal newStdGw=idLiftStd(newGw,&L);
164
165//what it looked like before idStd and idLift were replaced
166//by idLiftStd
167//   optionState=test;
168//   test|=Sy_bit(OPT_REDSB);
169//   test|=Sy_bit(OPT_NOT_BUCKETS);
170
171//   //PrintS("  new initial forms:\n");
172//   for (int ii=0; ii <IDELEMS(newGw); ii++) pCleardenom(newGw->m[ii]);
173//   //idShow(newGw);
174//   ideal newStdGw=idStd(newGw);
175//     PrintS("  std for initial forms done\n");
176
177//   test=optionState;
178
179//  matrix L=matIdLift(newGw,newStdGw);
180
181//  idDelete(&newStdGw);
182
183  idDelete(&newGw);
184    //PrintS("  lift for initial forms done\n");
185
186  ideal nextG=idrMoveR(G,oldRing,currRing);
187  rKill(oldRing);
188
189  matrix nextGmat=(matrix)nextG;
190
191  matrix resMat=mp_Mult(nextGmat,L,currRing);
192  idDelete((ideal *)&nextGmat);
193  idDelete((ideal *)&L);
194    //PrintS("  lift done\n");
195
196  nextG=(ideal)resMat;
197
198  optionState=test;
199  test|=Sy_bit(OPT_REDSB);
200  nextG = idInterRed(nextG);
201  test=optionState;
202
203  G=nextG;
204  return(state);
205}
206
207///////////////////////////////////////////////////////////////////
208
209
210///////////////////////////////////////////////////////////////////
211//walk64
212///////////////////////////////////////////////////////////////////
213//Description: the main function of groebner walk, keeping
214//everything together
215///////////////////////////////////////////////////////////////////
216//Uses: mStd, getnthrow, firstwalkStep64, nextt64,
217//nextw64, walkStep64, changecurrRing
218///////////////////////////////////////////////////////////////////
219
220WalkState walk64(ideal I,int64vec* currw64,ring destRing,
221int64vec* destVec64,ideal  & destIdeal,BOOLEAN sourceIsSB){
222
223  //some initializations
224  WalkState state=WalkOk;
225  BITSET optionState;
226  test|=Sy_bit(OPT_REDTAIL);
227  overflow_error=FALSE;
228  int step=0;
229  ideal G=I;
230
231  optionState=test;
232  test|=Sy_bit(OPT_REDSB);
233  if(!sourceIsSB)
234  {
235    ideal GG=idStd(G);
236    idDelete(&G); G=GG;
237  }
238  else
239    G=idInterRed(G);
240  test=optionState;
241
242  ideal nextG;
243  state=firstWalkStep64(G,currw64,destRing);
244  nextG=G;
245
246  if(overflow_error){
247      state=WalkOverFlowError;
248    return(state);
249    }
250
251   int64 nexttvec0,nexttvec1;
252   //int64vec* nexttvec64=nextt64(nextG,currw64,destVec64);
253   nextt64(nextG,currw64,destVec64,nexttvec0,nexttvec1);
254
255  //while(0<t<=1) ( t=((*nexttvec64)[0])/((*nexttvec64)[1]) )
256  //while( (*nexttvec64)[0]<=(*nexttvec64)[1] ) {
257  while (nexttvec0<=nexttvec1 ) {
258
259    step=step+1;
260
261    //int64vec *tt=nextw64(currw64,destVec64,nexttvec64);
262    int64vec *tt=nextw64(currw64,destVec64,nexttvec0,nexttvec1);
263    delete currw64; currw64=tt; tt=NULL;
264
265    if (TEST_OPT_PROT)
266    {
267      PrintS("walk step:"); currw64->show(); PrintLn();
268    }
269
270    state=walkStep64(nextG,currw64,step);
271    //uppdates nextG if all is OK
272
273    if(overflow_error){
274      return(WalkOverFlowError);
275    }
276
277    //delete nexttvec64;
278    //nexttvec64=nextt64(nextG,currw64,destVec64);
279    nextt64(nextG,currw64,destVec64,nexttvec0,nexttvec1);
280
281  }
282
283  destIdeal=sortRedSB(nextG);
284  return(state);
285}
286
287///////////////////////////////////////////////////////////////////
288
289
290
291///////////////////////////////////////////////////////////////////
292//FRACTAL WALK PART BEGINS HERE
293///////////////////////////////////////////////////////////////////
294
295///////////////////////////////////////////////////////////////////
296//firstFractalWalkStep64
297///////////////////////////////////////////////////////////////////
298//Description: called once before fractalRec64 is called, to set up
299//the ring and move G into it, has two different strategies
300///////////////////////////////////////////////////////////////////
301//Uses: firstWalkStep64,currwOnBorder64,getTaun64,rCopy0AndAddA,
302//getiv64,iv64Size,rComplete,rChangeCurrRing,idrMoveR,idStd,matIdLift
303///////////////////////////////////////////////////////////////////
304
305
306//unperturbedStartVectorStrategy IS NOW NOT ALLWAYS AS DEFAULT SET
307//TO TRUE BUT IS INPUT FROM fractalWalk64
308WalkState firstFractalWalkStep64(ideal & G,int64vec* & currw64,
309intvec* currMat, ring destRing,
310BOOLEAN unperturbedStartVectorStrategy){
311
312    //This strategy Uses the ordinary walk for the first step
313    if(unperturbedStartVectorStrategy){
314      return(unperturbedFirstStep64(G,currw64,destRing));
315    //here G is updated since its adress is given as argument
316    }
317
318    //This strategy makes sure that the start vector lies inside the start cone
319    //thus no step is needed within the start cone.
320    else{
321      if(currwOnBorder64(G,currw64))
322      {
323        int64 dummy64;
324        getTaun64(G,currMat,iv64Size(currw64),&currw64,dummy64);
325        //currw64=getiv64(getTaun64(G,currMat,iv64Size(currw64)));
326      }
327      ring oldRing=currRing;
328      ring newRing=rCopy0AndAddA(destRing,currw64);
329      rComplete(newRing);
330      rChangeCurrRing(newRing);
331      G=idrMoveR(G,oldRing,newRing);
332    }
333
334  return(WalkOk);
335}
336
337///////////////////////////////////////////////////////////////////
338
339
340
341//DIESE FUNKTION ERSETZT firstWalkStep FUR fractalWalk
342///////////////////////////////////////////////////////////////////
343//unperturbedFirstStep64
344///////////////////////////////////////////////////////////////////
345//Description: adapts currRing for the algorithm and moves G into
346//it, calculating the appropriate GB if currw not is inside of
347//the current groebner cone
348///////////////////////////////////////////////////////////////////
349//Assumes: that the option redSB is turned off
350///////////////////////////////////////////////////////////////////
351//Uses: init64, rCopy0AndAddA, rComplete, rChangeCurrRing, matIdLift,
352//mpMult, idInterRed, idStd, idCopy, idrMoveR,currwOnBorder64
353///////////////////////////////////////////////////////////////////
354
355WalkState unperturbedFirstStep64(ideal & G,int64vec* currw64, ring destRing){
356  WalkState state=WalkOk;
357  /* OLDRING **************************************************** */
358  ideal nextG;
359  BITSET optionState;
360
361  if (currwOnBorder64(G,currw64))
362  {
363    ideal Gw=init64(G,currw64);
364    ring oldRing=currRing;
365  /* NEWRING **************************************************** */
366    ring rnew=rCopy0AndAddA(destRing,currw64);
367    rComplete(rnew);
368    rChangeCurrRing(rnew);
369
370    ideal newGw=idrMoveR(Gw, oldRing,rnew);
371
372      //turn off bucket representation of polynomials and on redSB
373    optionState=test;
374    //test|=Sy_bit(OPT_NOT_BUCKETS);
375    test|=Sy_bit(OPT_REDSB);
376
377    ideal newStdGw=idStd(newGw);
378
379   //turn on bucket representation of polynomials and off redSB
380   test=optionState;
381
382    matrix L=matIdLift(newGw,newStdGw);
383    idDelete(&newStdGw);
384    idDelete(&newGw);
385
386    nextG=idrMoveR(G,oldRing,rnew); idTest(nextG);
387
388    matrix nextGmat=(matrix)nextG;
389
390    matrix resMat=mp_Mult(nextGmat,L,rnew);
391    idDelete((ideal *)&nextGmat);
392    idDelete((ideal *)&L);
393
394    nextG=(ideal)resMat;
395
396    optionState=test;
397    test|=Sy_bit(OPT_REDSB);
398    nextG = idInterRed(nextG);
399    test=optionState;
400  }
401  else
402  {
403    ring oldRing=currRing;
404    ring rnew=rCopy0AndAddA(destRing,currw64);
405    rComplete(rnew);
406    rChangeCurrRing(rnew);
407    nextG=idrMoveR(G,oldRing,rnew);
408  }
409
410  G=nextG;
411  return(state);
412}
413
414///////////////////////////////////////////////////////////////////
415
416
417
418///////////////////////////////////////////////////////////////////
419//fractalRec64
420///////////////////////////////////////////////////////////////////
421//Description: the recursion function used in fractalWalk64
422//(see appropriate par of thesis for more details)
423///////////////////////////////////////////////////////////////////
424//Assumes: that the option redSB is turned off
425///////////////////////////////////////////////////////////////////
426//Uses: getTaun64,getint64,getiv64,invEpsOk64,nextt64,nextw64,
427//init64,idCopy,noPolysOfMoreThanTwoTerms,rCopy0,rComplete,
428//rChangeCurrRing,rSetWeightVec,idrMoveR,mpMult,idInterRed
429///////////////////////////////////////////////////////////////////
430
431WalkState fractalRec64(ideal  & G,int64vec* currw64, intvec* destMat,
432                       int level,int step)
433{
434
435if (TEST_OPT_PROT)
436{  PrintS("fractal walk, weights");currw64->show();PrintLn(); }
437WalkState state=WalkOk;
438BITSET optionState=test;
439
440//1
441int64vec* w=(currw64);
442int64vec* old_w=(currw64);
443int64vec* sigma=(currw64);
444
445//lists taunpair64=getTaun64(G,destMat,level);
446//int64vec* w2=getiv64(taunpair64);
447//int64 inveps64=getint64(taunpair64);
448int64vec* w2;
449int64 inveps64;
450getTaun64(G,destMat,level,&w2,inveps64);
451
452//2
453while(1){
454
455  //int64vec* tvec64=nextt64(G,w,w2);
456  int64 tvec0,tvec1;
457  nextt64(G,w,w2,tvec0,tvec1);
458
459  if(overflow_error){
460    return WalkOverFlowError;
461  }
462
463  //tvec[1]>tvec[2] iff t>1 or t ist undefined i.e.is (2,0)
464  //if ((*tvec64)[0]>(*tvec64)[1])
465  if (tvec0>tvec1)
466  {
467    if(invEpsOk64(G,destMat,level,inveps64)) {
468      return(state);
469  }
470  else
471  {
472      //taunpair64=getTaun64(G,destMat,level);
473      //w2=getiv64(taunpair64);
474      //inveps64=getint64(taunpair64);
475      delete w2;
476      getTaun64(G,destMat,level,&w2,inveps64);
477
478      //tvec64=nextt64(G,w,w2);
479      nextt64(G,w,w2,tvec0,tvec1);
480
481      if(overflow_error){
482        return WalkOverFlowError;
483      }
484
485      //if((*tvec64)[0]>(*tvec64)[1])
486      if(tvec0>tvec1)
487      {
488        return(state);
489      }
490    }
491  }
492
493  //i.e. t=1 and we have reached the target vector
494  //if( ((*tvec64)[0]==(*tvec64)[1]) && (level!=iv64Size(w)) )
495  if( (tvec0==tvec1) && (level!=iv64Size(w)) )
496  {
497    state=fractalRec64(G,old_w,destMat,level+1,step);
498    return(state);
499  }
500  else
501  {
502    w=nextw64(w,w2,tvec0,tvec1);
503
504//3
505    ideal Gw=init64(G,w);  //finding initial ideal for new w
506
507    ring oldRing=currRing;
508
509    ideal GwCp=idCopy(Gw);
510    ideal GCp=idCopy(G);
511
512    ideal newGw;
513    ideal newStdGw;
514    ideal newG;
515
516//4
517    if ( level==iv64Size(w) || noPolysWithMoreThanTwoTerms(Gw) ){
518
519//5
520/*NEWRING**********************************************************/
521
522      //this assumes order has A-vector as first vector
523      ring newring=rCopy0(currRing);
524      rComplete(newring);
525      rSetWeightVec(newring,w->iv64GetVec());
526      rChangeCurrRing(newring);
527      //rSetWeightVec(newring,w->iv64GetVec());
528      //rComplete(newring,1);
529
530      newGw=idrMoveR(GwCp,oldRing,newring);
531
532      test|=Sy_bit(OPT_REDSB);
533      newStdGw=idStd(newGw); //computes new reduced GB of Gw
534      test=optionState;
535    }
536    else
537    {
538//7
539/*THE RING FROM THE RECURSION STEP BELOW***************************/
540      //Here we can choose whether to call fractalrec with old_w,
541      //which is the last w from the top level, or with sigma,
542      //the original start vector. The impact on the algorithm is not obvious.
543
544      state=fractalRec64(Gw,sigma,destMat,level+1,step);
545
546      //The resulting GB is Gw since its adress is given as argument.
547      ideal recG=Gw;
548      ring temp=currRing;
549
550/*NEWRING**********************************************************/
551
552      //this assumes order has A-vector as first vector
553      ring newring=rCopy0(currRing);
554      rComplete(newring);
555      rChangeCurrRing(newring);
556      rSetWeightVec(currRing,w->iv64GetVec());
557      rComplete(newring,1);
558
559      newGw=idrMoveR(GwCp,oldRing,newring);
560
561      newStdGw=idrMoveR(recG,temp,newring);
562    }
563
564//8
565    //lifting comes either after level=nvars(ring), after Gw has
566    //no poly with more than two terms or after
567    //fractalRec64(Gw,level+1) has returned
568
569    //test|=Sy_bit(OPT_NOT_BUCKETS);
570    matrix L=matIdLift(newGw,newStdGw);
571    test=optionState;
572
573    newG=idrMoveR(GCp,oldRing,currRing);
574    matrix MG=(matrix)newG;
575
576    matrix resMat=mp_Mult(MG,L,currRing);
577    idDelete((ideal *)&MG);
578    idDelete((ideal *)&L);
579    G=(ideal)resMat;
580
581//9
582
583    test|=Sy_bit(OPT_REDSB);
584    G=idInterRed(G);
585    test=optionState;
586
587    old_w=iv64Copy(w);
588    if(level==1) step=step+1;
589
590  }
591
592}
593}
594
595
596///////////////////////////////////////////////////////////////////
597
598
599//ANOTHER INPUT-VARIABLE ADDED: unperturbedStartVectorStrategy
600//THIS SHOULD BE SET IN walkProc.cc BY THE USER
601///////////////////////////////////////////////////////////////////
602//fractalWalk64
603///////////////////////////////////////////////////////////////////
604//Description:
605///////////////////////////////////////////////////////////////////
606//Uses:fractalRec64,firstFractalWalkStep64,idStd,idInterRed,
607//int64VecToIntVec,getNthRow,sortRedSB
608///////////////////////////////////////////////////////////////////
609
610WalkState fractalWalk64(ideal sourceIdeal,ring destRing,
611ideal & destIdeal,BOOLEAN sourceIsSB,
612BOOLEAN unperturbedStartVectorStrategy)
613{
614
615  overflow_error=FALSE; //global
616  WalkState state=WalkOk;
617  test|=Sy_bit(OPT_REDTAIL);
618  ideal G;
619
620  BITSET optionState=test;
621  test|=Sy_bit(OPT_REDSB);
622  if(!sourceIsSB)
623  {
624    G=idStd(sourceIdeal);
625  }
626
627  else
628  {
629    G=idInterRed(idCopy(sourceIdeal));
630  }
631
632  test=optionState; //switches REDSB off
633
634  //matrices for the orders of the rings
635  intvec *destMat=int64VecToIntVec(rGetGlobalOrderMatrix(destRing));
636  intvec *currMat=int64VecToIntVec(rGetGlobalOrderMatrix(currRing));
637
638  int64vec* currw64=getNthRow64(currMat,1); //start vector
639
640  state=firstFractalWalkStep64(G,currw64,currMat,destRing,
641                               unperturbedStartVectorStrategy);
642  delete currMat;
643
644  state=fractalRec64(G,currw64,destMat,1,1);
645  if(state==WalkOk)
646    destIdeal=G;
647
648  if(overflow_error)
649    state=WalkOverFlowError;
650
651  delete currw64;
652  delete destMat;
653  return state;
654}
655
656
657/////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.