source: git/kernel/walkMain.cc @ fbc7cb

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