Changeset 5684df in git for kernel/gfan.cc

Timestamp:

May 4, 2009, 5:14:38 PM (15 years ago)

Author:

Martin Monerjan

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

b2a6e3e0d28bf571be753914eaccf9f5c629080f

Parents:

e2efbe929b401e118eadf323f44b5bd894b457b0

Message:

found bug in flippability test


git-svn-id: file:///usr/local/Singular/svn/trunk@11767 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• kernel/gfan.cc (modified) (10 diffs)

kernel/gfan.cc

@@ -2 +2 @@
 Compute the Groebner fan of an ideal
 $Author: monerjan $
-$Date: 2009-04-30 09:32:01 $
-$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.43 2009-04-30 09:32:01 monerjan Exp $
-$Id: gfan.cc,v 1.43 2009-04-30 09:32:01 monerjan Exp $
+$Date: 2009-05-04 15:14:38 $
+$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.44 2009-05-04 15:14:38 monerjan Exp $
+$Id: gfan.cc,v 1.44 2009-05-04 15:14:38 monerjan Exp $
 */
 
@@ -210 +210 @@
                 
                 /** \brief Default destructor */
-                ~gcone();               //destructor
+                ~gcone(){;}             //destructor
                 
                 /** Pointer to the first facet */
@@ -403 +403 @@
                                 }//for (int jj = 1; jj <ddcols; jj++)
                                 /*Quick'n'dirty hack for flippability*/ 
-                                bool isFlippable;                       
+                                bool isFlippable;
+                                //NOTE BUG HERE
+                                /* The flippability check is wrong!
+                                (1,-4) will pass, but (-1,7) will not. 
+                                REWRITE THIS
+                                */                      
                                 for (int jj = 0; jj<this->numVars; jj++)
                                 {                                       
@@ -514 +519 @@
                                         cout << endl;
 #endif
-                                        if (isParallel(check,fNormal)) //pass *check when 
+                                        //TODO why not *check, *fNormal????
+                                        if (isParallel(*check,*fNormal)) //pass *check when 
                                         {
                                                 cout << "Parallel vector found, adding to initialFormElement" << endl;                  
@@ -1022 +1028 @@
                 /** 
                 * Determines whether a given facet of a cone is the search facet of a neighbouring cone
-                *
+                * This is done in the following way:
+                * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
+                * that is first crossed during the generic walk.
+                * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
+                * If this is the case, then our facet is indeed a search facet and TRUE is retuned. 
                 */
-                bool isSearchFacet(gcone &tmpcone, facet *testfacet)
-                {                       
+                bool isSearchFacet(gcone &tmpcone, facet &testfacet)
+                {                               
                         ring actRing=currRing;
-                        
-                        facet *fMin=new facet();        //Pointer to the "minimal" facet
+                        facet *facetPtr=(facet*)tmpcone.facetPtr;                       
+                        facet *fMin=new facet(*facetPtr);       //Pointer to the "minimal" facet
+                        //facet *fMin = new facet(tmpcone.facetPtr);
+                        //fMin=tmpcone.facetPtr;                //Initialise to first facet of tmpcone
                         facet *fAct;
-                        fMin = testfacet;
+                        //fMin = testfacet;
+                        //fMin->next=testfacet->next;
                         fAct = fMin;
                         
-                        rChangeCurrRing(this->rootRing);
+                        cout << endl<< fMin->next << endl;
+                        cout << tmpcone.facetPtr->next << endl;
+                        tmpcone.facetPtr->printNormal();
+                        fAct=tmpcone.facetPtr->next;
+                        fAct->printNormal();
+                        /*fMin->printNormal();
+                        fMin->next->printNormal();*/
+                        
+                        rChangeCurrRing(this->rootRing);        //because we compare the monomials in the rootring                      
                         poly p=NULL;
                         poly q=NULL;
@@ -1040 +1061 @@
                         intvec *sigma = new intvec(this->numVars);
                         sigma=tmpcone.getIntPoint();
+                        cout << "pling" << endl;
                         int *u=(int *)omAlloc((this->numVars+1)*sizeof(int));
                         int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
                         int weight1,weight2;
-                        while(tmpcone.facetPtr->next!=NULL)
+                        while(fAct->next!=NULL) //ersetzen durch fAct
                         {
                                 /* Get alpha_i and alpha_{i+1} */
                                 p_weight=fMin->getFacetNormal();
                                 q_weight=fMin->next->getFacetNormal();
-                                
+                                p_weight->show(); 
+                                q_weight->show();
                                 /*Compute the dot product of sigma and alpha_{i,j}*/
                                 weight1=dotProduct(sigma,p_weight);
                                 weight2=dotProduct(sigma,q_weight);
-                                
+                                cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
                                 /*Adjust alpha_i and alpha_i+1 accordingly*/
                                 for(int ii=1;ii<=this->numVars;ii++)
@@ -1061 +1084 @@
                                         v[ii]=weight2*(*q_weight)[ii];
                                 }
+                                cout << "PLONK" << endl;
                                 
                                 /*Now p_weight and q_weight need to be compared as exponent vectors*/
                                 pSetExpV(p,u);
                                 pSetExpV(q,v);
+                                cout << "AARGH" << endl;
                                 /*We want to check whether x^p < x^q 
                                 => want to check for return value 1 */
@@ -1076 +1101 @@
                                         fAct=fAct->next;
                                 }
+                                fAct=fAct->next;
                         }//while(tmpcone.facetPtr->next!=NULL)
-                        rChangeCurrRing(actRing);
+                        
+                        /*If testfacet was minimal then fMin should still point there */
+                        //if (isParallel(fMin->getFacetNormal,testfacet.getFacetNormal))
+                        if (fMin==tmpcone.facetPtr)
+                        {
+                                rChangeCurrRing(actRing);
+                                return TRUE;
+                        }
+                        else 
+                        {
+                                rChangeCurrRing(actRing);
+                                return FALSE;
+                        }
                 }//bool isSearchFacet
                 
-                void reverseSearch(gcone gcAct) //no const possible here since we call gcAct->flip
+                void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
                 {
                         facet *fAct=new facet();
-                        fAct = gcAct.facetPtr;                  
-                        
-                        while(fAct->next!=NULL)
-                        {
-                                gcAct.flip(gcAct.gcBasis,gcAct.facetPtr);
-                                gcone *gcTmp = new gcone(gcAct);
+                        fAct = gcAct->facetPtr;                 
+                        
+                        while(fAct->next!=NULL)  //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
+                        {
+                                cout << "==========================================================================================="<< endl;
+                                gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
+                                gcone *gcTmp = new gcone(*gcAct);
                                 idShow(gcTmp->gcBasis);
                                 gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
+#ifdef gfan_DEBUG
+                                facet *f = new facet();
+                                f=gcTmp->facetPtr;
+                                while(f->next!=NULL)
+                                {
+                                        f->printNormal();
+                                        f=f->next;                                      
+                                }
+#endif
                                 gcTmp->showIntPoint();
                                 /*recursive part goes gere*/
-                                if (isSearchFacet(gcTmp,gcAct.facetPtr))
-                                {
-                                        gcAct.next=gcTmp;
-                                        reverseSearch(*gcTmp);
+                                if (isSearchFacet(*gcTmp,(facet&)gcAct->facetPtr))
+                                {
+                                        gcAct->next=gcTmp;
+                                        cout << "PING"<< endl;
+                                        reverseSearch(gcTmp);
                                 }
                                 else
@@ -1113 +1162 @@
         int numvar = pVariables; 
         
-        #ifdef gfan_DEBUG
+#ifdef gfan_DEBUG
         cout << "Now in subroutine gfan" << endl;
-        #endif
+#endif
         ring inputRing=currRing;        // The ring the user entered
         ring rootRing;                  // The ring associated to the target ordering
@@ -1179 +1228 @@
                 fAct = fAct->next;              
         }*/
-        gcAct->reverseSearch(*gcAct);
+        gcAct->reverseSearch(gcAct);
         /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
         The return type will then be of type LIST_CMD