Changeset f5b077 in git for kernel/gfan.cc
 Timestamp:
 Apr 6, 2009, 4:57:18 PM (14 years ago)
 Branches:
 (u'spielwiese', '91e5db82acc17434e4062bcfa44e6efa7d41fd30')
 Children:
 1288efb1bc6f003f059484c99a625cf47c1ed08b
 Parents:
 f3c31a15740169790d1f07f8cb43045dd6f97991
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

kernel/gfan.cc
rf3c31a rf5b077 2 2 Compute the Groebner fan of an ideal 3 3 $Author: monerjan $ 4 $Date: 2009040 3 13:49:16$5 $Header: /exports/cvsroot2/cvsroot/kernel/gfan.cc,v 1.2 8 20090403 13:49:16monerjan Exp $6 $Id: gfan.cc,v 1.2 8 20090403 13:49:16monerjan Exp $4 $Date: 20090406 14:57:18 $ 5 $Header: /exports/cvsroot2/cvsroot/kernel/gfan.cc,v 1.29 20090406 14:57:18 monerjan Exp $ 6 $Id: gfan.cc,v 1.29 20090406 14:57:18 monerjan Exp $ 7 7 */ 8 8 … … 52 52 { 53 53 private: 54 /** inner normal, describing the facet uniquely up to isomorphism */ 55 intvec *fNormal; 54 /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */ 55 intvec *fNormal; 56 57 /** \brief The Groebner basis on the other side of a shared facet 58 * 59 * In order not to have to compute the flipped GB twice we store the basis we already get 60 * when identifying search facets. Thus in the next step of the reverse search we can 61 * just copy the old cone and update the facet and the gcBasis. 62 * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB 63 */ 64 ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip 65 66 56 67 public: 57 68 /** The default constructor. Do I need a constructor of type facet(intvec)? */ … … 85 96 } 86 97 87 /** \brief The Groebner basis on the other side of a shared facet 88 * 89 * In order not to have to compute the flipped GB twice we store the basis we already get 90 * when identifying search facets. Thus in the next step of the reverse search we can 91 * just copy the old cone and update the facet and the gcBasis 92 */ 93 ideal flibGB; //The Groebner Basis on the other side, computed via gcone::flip 98 /** Store the flipped GB*/ 99 void setFlipGB(ideal I) 100 { 101 this>flipGB=I; 102 } 103 104 /** Print the flipped GB*/ 105 void printFlipGB() 106 { 107 idShow(this>flipGB); 108 } 94 109 95 110 bool isFlippable; //flippable facet? Want to have cone>isflippable.facet[i] … … 143 158 * As a result of this procedure the pointer facetPtr points to the first facet of the cone. 144 159 */ 145 void getConeNormals( idealI)160 void getConeNormals(const ideal &I) 146 161 { 147 162 #ifdef gfan_DEBUG … … 330 345 /*1st step: Compute the initial ideal*/ 331 346 poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form 332 ideal initialForm ;347 ideal initialForm=idInit(IDELEMS(gb),this>gcBasis>rank); 333 348 poly aktpoly; 334 349 intvec *check = new intvec(this>numVars); //array to store the difference of LE and v … … 340 355 int *leadExpV=(int *)omAlloc((this>numVars+1)*sizeof(int)); 341 356 pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p) 342 //pGetExpV(aktpoly,ivLeadExpV);343 357 initialFormElement[ii]=pHead(aktpoly); 344 358 … … 355 369 (*check)[jj]=v[jj+1]leadExpV[jj+1]; 356 370 } 371 #ifdef gfan_DEBUG 357 372 cout << "check="; 358 373 check>show(); 359 374 cout << endl; 375 #endif 360 376 if (isParallel(check,fNormal)) //pass *check when 361 377 { … … 367 383 pWrite(initialFormElement[ii]); 368 384 cout << "" << endl; 369 /*Now initialFormElement must be added to (ideal)initialForm */ 385 /*Now initialFormElement must be added to (ideal)initialForm */ 386 //f>flibGB>m[ii]=(poly)initialFormElement[ii]; 387 //(poly)initialForm>m[ii]=pAdd(initialForm[ii],initialFormElement[ii]); 388 initialForm>m[ii]=initialFormElement[ii]; 370 389 }//for 390 f>setFlipGB(initialForm); 391 #ifdef gfan_DEBUG 392 cout << "Initial ideal is: " << endl; 393 //idShow(initialForm); 394 f>printFlipGB(); 395 cout << "===" << endl; 396 #endif 397 delete check; 398 399 /*2nd step: lift initial ideal to a GB of the neighbouring cone*/ 400 ring tmpring=rCopy0(currRing); 401 tmpring>order[0]=ringorder_a; 402 tmpring>order[1]=ringorder_dp; 403 tmpring>order[2]=ringorder_C; 404 //rWrite(tmpring); 405 tmpring>wvhdl[0] =( int *)omAlloc((fNormal>length())*sizeof(int)); //found in Singular/ipshell.cc 406 407 for (int ii=0;ii<this>numVars;ii++) 408 { 409 cout << "ping" << endl; 410 tmpring>wvhdl[0][ii]=(*fNormal)[ii]; //What exactly am I doing here? 411 cout << "pong" << endl; 412 cout << *tmpring>wvhdl[ii] << endl; 413 } 414 //tmpring>wvhdl=(int**)(fNormal); 415 rComplete(tmpring); 416 rWrite(tmpring); 417 /*setring(tmpring); 418 ideal ina=initialForm; 419 ideal H; 420 H=kstd(ina,NULL,testHomog,NULL); 421 idShow(H);*/ 422 371 423 }//void flip(ideal gb, facet *f) 372 424 … … 377 429 *\return void 378 430 */ 379 void getGB( idealinputIdeal)431 void getGB(const ideal &inputIdeal) 380 432 { 381 433 ideal gb; … … 410 462 bool isParallel(intvec *a, intvec *b) 411 463 { 412 //bool res=FALSE;413 //intvec iva(this>numVars);414 //intvec ivb(this>numVars);415 464 int lhs,rhs; 416 //lhs=0;417 //rhs=0;418 //iva = a;419 //ivb = b;420 //lhs=dotProduct(&iva,&ivb)*dotProduct(&iva,&ivb);421 //lhs=dotProduct(&iva,&iva)*dotProduct(&ivb,&ivb);422 465 lhs=dotProduct(&a,&b)*dotProduct(&a,&b); 423 466 rhs=dotProduct(&a,&a)*dotProduct(&b,&b); … … 476 519 3. getConeNormals 477 520 */ 478 479 521 480 522 /* Construct a new ring which will serve as our root … … 482 524 */ 483 525 rootRing=rCopy0(currRing); 484 rootRing>order[0]=ringorder_ dp;526 rootRing>order[0]=ringorder_lp; 485 527 rComplete(rootRing); 486 //rChangeCurrRing(rootRing); 528 rChangeCurrRing(rootRing); 529 530 /* Fetch the inputIdeal into our rootRing */ 531 map theMap=(map)idMaxIdeal(1); //evil hack! 532 //idShow(idMaxIdeal(1)); 533 /*for (int ii=0;ii<pVariables;ii++) 534 { 535 theMap>m[ii]=inputIdeal>m[ii]; 536 }*/ 537 theMap>preimage=NULL; 487 538 ideal rootIdeal; 488 /* Fetch the inputIdeal into our rootRing */ 489 map m=(map)idInit(IDELEMS(inputIdeal),0); 490 //rootIdeal=fast_map(inputIdeal,inputRing,(ideal)m, currRing); 539 rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing); 491 540 #ifdef gfan_DEBUG 492 541 cout << "Root ideal is " << endl; 493 id Print(rootIdeal);494 cout << "The current ring is " << endl;542 idShow(rootIdeal); 543 cout << "The root ring is " << endl; 495 544 rWrite(rootRing); 496 545 cout << endl; … … 501 550 gcAct = gcRoot; 502 551 gcAct>numVars=pVariables; 503 gcAct>getGB(inputIdeal); 552 gcAct>getGB(rootIdeal); //sets gcone::gcBasis 553 idShow(gcAct>gcBasis); 504 554 gcAct>getConeNormals(gcAct>gcBasis); //hopefully compute the normals 505 555 gcAct>flip(gcAct>gcBasis,gcAct>facetPtr);
Note: See TracChangeset
for help on using the changeset viewer.