Changeset 73809b in git

Timestamp:

Mar 6, 2010, 2:36:13 PM (13 years ago)

Author:

Martin Monerjan

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

e6ddbdacf047de338a179b12ccb98096d4f02871

Parents:

c57a83edd49526195f04ddf602fca071cafed3b7

Message:

#defined SHALLOW
facet::shallowCopy, facet::shallowDelete. Needed for gcone::enqueue2
Fixes in gcone::areEqual2(facet*,facet*), now the standard for
comparisons of facet
facet::getRef2InteriorPoint
getExtremalRays: Fixed dividing out of gcd
makeInt: now returns primitve vector
renamed gcone::areEqual(intvec&, intvec&) to gcone::ivAreEqual
gcone::normalize no longer needed



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

Location:

kernel

Files:

• gfan.cc (modified) (59 diffs)
• gfan.h (modified) (8 diffs)

kernel/gfan.cc

@@ -75 +75 @@
 #endif
 
+//NOTE Defining this will slow things down!
+//Only good for very coarse profiling
 #define gfanp
 #ifdef gfanp
 #include <sys/time.h>
+#endif
+
+#ifndef SHALLOW
+#define SHALLOW
 #endif
 
@@ -114 +120 @@
 * could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
 */
-facet::facet(int const &n)
+facet::facet(const int &n)
 {
         this->fNormal=NULL;
@@ -144 +150 @@
         //Needed for flip2
         //NOTE ONLY REFERENCE
-        this->interiorPoint=/*ivCopy(*/f.interiorPoint/*)*/;//only referencing is prolly not the best thing to do in a copy constructor
+        this->interiorPoint=ivCopy(f.interiorPoint);//only referencing is prolly not the best thing to do in a copy constructor
         facet* f2Copy;
         f2Copy=f.codim2Ptr;
@@ -151 +157 @@
         while(f2Copy!=NULL)
         {
-                if(f2Act==NULL)
+                if(f2Act==NULL || f2Act==(facet*)0xfefefefefefefefe)
                 {
                         f2Act=new facet(2);
@@ -175 +181 @@
 
 /** \brief Shallow copy constructor for facets
+* We only need the interior point for equality testing
 */
-facet::facet(const facet& f, bool shallow)
-{}
+facet* facet::shallowCopy(const facet& f)
+{
+        facet *res = new facet();
+        res->fNormal=(intvec* const)f.fNormal;
+        res->UCN=f.UCN;
+        res->isFlippable=f.isFlippable;
+        res->interiorPoint=(intvec * const)f.interiorPoint;
+        res->codim2Ptr=(facet * const)f.codim2Ptr;
+        res->prev=NULL;
+        res->next=NULL;
+        res->flipGB=NULL;
+        return res;
+}
+
+void facet::shallowDelete()
+{
+        this->fNormal=NULL;
+        this->UCN=0;
+        this->interiorPoint=NULL;
+        this->codim2Ptr=NULL;
+        this->prev=NULL;
+        this->next=NULL;
+        this->flipGB=NULL;
+//      delete this;
+}
                 
 /** The default destructor */
@@ -217 +247 @@
         return(this->fNormal);
 }       
-        
-bool gcone::areEqual2(facet* f, facet *g)
+
+/** Equality check for facets based on unique interior points*/
+inline bool gcone::areEqual2(facet* f, facet *g)
 {
 #ifdef gfanp
@@ -225 +256 @@
         gettimeofday(&start, 0);
 #endif
-        bool res = FALSE;
-        const intvec* fNormal;
-        const intvec* gNormal;
-        fNormal = f->getRef2FacetNormal();
-        gNormal = g->getRef2FacetNormal();
-        intvec *fNRef=const_cast<intvec*>(fNormal);
-        intvec *gNRef=const_cast<intvec*>(gNormal);
-        if(isParallel(*fNRef,*gNRef))
-        {
+        bool res = TRUE;
+//      const intvec *fNormal;
+//      const intvec *gNormal;
+//      fNormal = f->getRef2FacetNormal();
+//      gNormal = g->getRef2FacetNormal();
+//      intvec *fNRef=const_cast<intvec*>(fNormal);
+//      intvec *gNRef=const_cast<intvec*>(gNormal);
+        /*if(isParallel(fNRef,*gNRef))
+        {
+                const intvec *fIntP = f->getRef2InteriorPoint();
+                const intvec *gIntP = g->getRef2InteriorPoint();
                 for(int ii=0;ii<this->numVars;ii++)
                 {
-                        if( (*fNormal)[ii]!=(*gNormal)[ii] )
-                        {
-                                res=TRUE;
+                        if( (*fIntP)[ii] != (*gIntP)[ii] )
+                        {
+                                res=FALSE;
                                 break;
                         }
                 }
         }
-        return res;
+        else
+                res=FALSE;*/
+        const intvec *fIntP = f->getRef2InteriorPoint();
+        const intvec *gIntP = g->getRef2InteriorPoint();
+        for(int ii=0;ii<this->numVars;ii++)
+        {
+                if( (*fIntP)[ii] != (*gIntP)[ii] )
+                {
+                        res=FALSE;
+                        break;
+                }
+        }
 #ifdef gfanp
         gettimeofday(&end, 0);
         t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
 #endif  
+        return res;
 }
 
@@ -385 +430 @@
 //      return this->fNormal;
 }
-                
+
 /** Method to print the facet normal*/
 inline void facet::printNormal()
@@ -456 +501 @@
 {
         return ivCopy(this->interiorPoint);
-}       
+}
+
+inline const intvec *facet::getRef2InteriorPoint()
+{
+        return (this->interiorPoint);
+}
                 
 /** \brief Debugging function
@@ -565 +615 @@
 gcone::~gcone()
 {
+#ifndef NDEBUG
+  #if SIZEOF_LONG==8
+        if(this->gcBasis!=(ideal)(0xfbfbfbfbfbfbfbfb))
+                idDelete((ideal*)&this->gcBasis);
+  #elif SIZEOF_LONG!=8
+        if(this->gcBasis!=(ideal)0xfbfbfbfb)
+                idDelete((ideal *)&this->gcBasis);
+  #endif
+#else
         if(this->gcBasis!=NULL)
                 idDelete((ideal *)&this->gcBasis);
+#endif
+//      idDelete((ideal *)&this->gcBasis);
 //      if(this->inputIdeal!=NULL)
 //              idDelete((ideal *)&this->inputIdeal);
@@ -589 +650 @@
         //dd_FreeMatrix(this->ddFacets);
 }                       
-        
+
+/** Returns the number of cones existing at the time*/
 inline int gcone::getCounter()
 {
@@ -603 +665 @@
 }
                 
-/** \brief Return the interior point */
-inline intvec *gcone::getIntPoint()
-{
-        return ivCopy(this->ivIntPt);
+/** \brief Returns either a physical copy the interior point of a cone or just a reference to it.*/
+inline intvec *gcone::getIntPoint(bool shallow)
+{
+        if(shallow==TRUE)
+                return this->ivIntPt;
+        else
+                return ivCopy(this->ivIntPt);
 }
                 
@@ -762 +827 @@
 {
         return this->pred;
+}
+/** Returns a copy of the this->baseRing */
+inline ring gcone::getBaseRing()
+{
+        return rCopy(this->baseRing);
 }
 
@@ -796 +866 @@
         for (int ii=0;ii<IDELEMS(I);ii++)
         {
-                aktpoly=(poly)I->m[ii];
-                rows=rows+pLength(aktpoly)-1;
+//              aktpoly=(poly)I->m[ii];
+//              rows=rows+pLength(aktpoly)-1;
+                rows=rows+pLength((poly)I->m[ii])-1;
         }
 
@@ -1005 +1076 @@
         for (int kk = 0; kk<ddrows; kk++)
         {
+                int ggT=1;//NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
                 intvec *load = new intvec(this->numVars);//intvec to store a single facet normal that will then be stored via setFacetNormal
                 for (int jj = 1; jj <ddcols; jj++)
@@ -1010 +1082 @@
                         double foo;
                         foo = mpq_get_d(ddineq->matrix[kk][jj]);
-                        (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load            
+                        (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load
+                        ggT = intgcd(ggT,(int&)foo);
                 }//for (int jj = 1; jj <ddcols; jj++)
-                                
+                if(ggT>1)
+                {
+                        for(int ll=0;ll<this->numVars;ll++)
+                                (*load)[ll] /= ggT;//make primitive vector
+                }
                 /*Quick'n'dirty hack for flippability. Executed only if gcone::hasHomInput==FALSE
                 * Otherwise every facet intersects the positive orthant
@@ -1312 +1389 @@
 * So we keep pointer tmp to iv and delete(tmp), so there should not occur a 
 * memleak
+* TODO normalization
 */
 void gcone::getExtremalRays(const gcone &gc)
@@ -1318 +1396 @@
         timeval start, end;
         gettimeofday(&start, 0);
+        timeval poly_start, poly_end;
+        gettimeofday(&poly_start,0);
 #endif
         //Add lineality space - dd_LinealitySpace
@@ -1347 +1427 @@
         dd_MatrixPtr P;
         P=dd_CopyGenerators(ddPolyh);
-        /*dd_rowset ddlinset,ddredrows; dd_rowindex ddnewpos;
-        dd_MatrixCanonicalize(&P, &ddlinset, &ddredrows, &ddnewpos, &err);*/
         dd_FreePolyhedra(ddPolyh);
+        dd_FreeMatrix(ddineq);
+#ifdef gfanp
+        gettimeofday(&poly_end,0);
+        t_ddPolyh += (poly_end.tv_sec - poly_start.tv_sec + 1e-6*(poly_end.tv_usec - poly_start.tv_usec));
+#endif
         /* Compute interior point on the fly*/
         intvec *ivIntPointOfCone = new intvec(this->numVars);
+//      for(int ii=0;ii<P->rowsize;ii++)
+//      {}
         mpq_t *colSum = new mpq_t[this->numVars];
         int denom[this->numVars];//denominators of colSum
+        //NOTE TODO need to gcd of rows and factor out! -> makeInt
         for(int jj=0;jj<this->numVars;jj++)
         {
@@ -1372 +1458 @@
         }
         //Now compute lcm of denominators of colSum[jj]
+        //NOTE gcd as well and normalise instantly?
         mpz_t kgV; mpz_init(kgV);
         mpz_set_str(kgV,"1",10);
@@ -1389 +1476 @@
         mpz_clear(den);
         mpz_clear(tmp);
+        int ggT=(*ivIntPointOfCone)[0];
         for (int ii=0;ii<(this->numVars);ii++)
         {
@@ -1396 +1484 @@
                 (*ivIntPointOfCone)[ii]=(int)mpz_get_d(mpq_numref(product));
                 mpq_clear(product);
+                //Compute intgcd
+                ggT=intgcd(ggT,(*ivIntPointOfCone)[ii]);
+        }
+        //Divide out a common gcd > 1
+        if(ggT>1)
+        {
+                for(int ii=0;ii<this->numVars;ii++)
+                        (*ivIntPointOfCone)[ii] /= ggT;
         }
         mpq_clear(qkgV);
+        delete [] colSum;
         /*For homogeneous input (like Det3,3,5) the int points may be negative. So add a suitable multiple of (1,_,1)*/
         if(hasHomInput==TRUE && iv64isStrictlyPositive(ivIntPointOfCone)==FALSE)
@@ -1413 +1510 @@
                 }
                 delete ivOne;
+                int ggT=(*ivIntPointOfCone)[0];
+                for(int ii=0;ii<this->numVars;ii++)
+                        ggT=intgcd( ggT, (*ivIntPointOfCone)[ii]);
+                if(ggT>1)
+                {
+                        for(int jj=0;jj<this->numVars;jj++)
+                                (*ivIntPointOfCone)[jj] /= ggT;
+                }
         }
         assert(iv64isStrictlyPositive(ivIntPointOfCone));
+        
         this->setIntPoint(ivIntPointOfCone);
         delete(ivIntPointOfCone);
@@ -1430 +1536 @@
                 {                       
                         intvec *rowvec = new intvec(this->numVars);
-                        makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational                      
+                        makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve                  
                         if(dotProduct(*fNormal,*rowvec)==0)
                         {
@@ -1456 +1562 @@
                         }
                         delete(rowvec);
-                }
+                }//For non-homog input ivIntPointOfFacet should already be >0 here
+                if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
                 //if we have no strictly pos ray but the input is homogeneous
                 //then add a suitable multiple of (1,...,1)
@@ -1477 +1584 @@
                                 delete tmp;                             
                         }
-                        fAct->setInteriorPoint(ivIntPointOfFacet);
+                        //fAct->setInteriorPoint(ivIntPointOfFacet);
                         delete ivOne;
                 }
-                else
-                        fAct->setInteriorPoint(ivIntPointOfFacet);
+                //else
+                int ggT=(*ivIntPointOfFacet)[0];
+                for(int ii=0;ii<this->numVars;ii++)
+                        ggT=intgcd(ggT,(*ivIntPointOfFacet)[ii]);
+                if(ggT>1)
+                {
+                        for(int ii=0;ii<this->numVars;ii++)
+                                (*ivIntPointOfFacet)[ii] /= ggT;
+                }                       
+                fAct->setInteriorPoint(ivIntPointOfFacet);
                 
                 delete(ivIntPointOfFacet);
@@ -1637 +1752 @@
         if( (srcRing->order[0]!=ringorder_a))
         {
-                intvec *iv = new intvec(this->numVars);
-                iv = ivNeg(fNormal);
+                intvec *iv;// = new intvec(this->numVars);
+                iv = ivNeg(fNormal);//ivNeg uses ivCopy -> new
 //              tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
                 tmpRing=rCopyAndAddWeight(srcRing,iv);
@@ -2345 +2460 @@
 inline int gcone::dotProduct(const intvec &iva, const intvec &ivb)                              
 {                       
-        int res=0;
+        int res=0;      
         for (int i=0;i<this->numVars;i++)
         {
@@ -2360 +2475 @@
  */
 inline bool gcone::isParallel(const intvec &a,const intvec &b)
-{                       
+{       
+/*#ifdef gfanp  
+        timeval start, end;
+        gettimeofday(&start, 0);
+#endif*/                
         int lhs,rhs;
         bool res;
@@ -2374 +2493 @@
                 res = FALSE;
         }
+// #ifdef gfanp
+//      gettimeofday(&end, 0);
+//      t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
+// #endif       
         return res;
 }//bool isParallel
@@ -2952 +3075 @@
 }
                 
+//NOTE not needed anywhere
 ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
 {
@@ -2973 +3097 @@
 /** \brief Check for equality of two intvecs
  */
-inline bool gcone::areEqual(intvec &a, intvec &b)
+inline bool gcone::ivAreEqual(intvec &a, intvec &b)
 {
         bool res=TRUE;
@@ -3102 +3226 @@
 }       
 
-/** Compute the kernel of a given mxn-matrix */
-// dd_MatrixPtr gcone::kernel(const dd_MatrixPtr &M)
-// {
-// }
 
 /** \brief The new method of Markwig and Jensen
@@ -3154 +3274 @@
                                 
         //Compute unique representation of Facets and rays, i.e. primitive vectors
-        gcAct->normalize();
+//      gcAct->normalize();
                         
         /* Make a copy of the facet list of first cone
@@ -3167 +3287 @@
                 if(fAct->isFlippable==TRUE)
                 {
+                        /*Using shallow copy here*/
+#ifdef SHALLOW
+                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
+                        {
+                                if(SearchListRoot!=NULL)
+                                        delete(SearchListRoot);
+                                SearchListRoot = fAct->shallowCopy(*fAct);
+                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
+                        }
+                        else
+                        {                       
+                                SearchListAct->next = fAct->shallowCopy(*fAct);
+                                SearchListAct = SearchListAct->next;                                            
+                        }
+                        fAct=fAct->next;
+#else
+                        /*End of shallow copy*/
                         intvec *fNormal;
                         fNormal = fAct->getFacetNormal();
@@ -3217 +3354 @@
                         }
                         fAct = fAct->next;
-                        delete fNormal;                 
+                        delete fNormal;
+#endif
                 }//if(fAct->isFlippable==TRUE)                  
                 else {fAct = fAct->next;}
@@ -3224 +3362 @@
         SearchListAct = SearchListRoot; //Set to beginning of list
         /*Make SearchList doubly linked*/
+#ifndef NDEBUG
+  #if SIZEOF_LONG==8
+        while(SearchListAct!=(facet*)0xfefefefefefefefe && SearchListAct!=NULL)
+        {
+                if(SearchListAct->next!=(facet*)0xfefefefefefefefe && SearchListAct->next!=NULL){
+  #elif SIZEOF_LONG!=8
+        while(SearchListAct!=(facet*)0xfefefefe)
+        {
+                if(SearchListAct->next!=(facet*)0xfefefefe){
+  #endif
+#else
         while(SearchListAct!=NULL)
         {
-                if(SearchListAct->next!=NULL)
-                {
+                if(SearchListAct->next!=NULL){
+#endif
+//              if(SearchListAct->next!=NULL)
+//              {
                         SearchListAct->next->prev = SearchListAct;                                      
                 }
@@ -3243 +3394 @@
           We always choose the first facet from SearchList as facet to be flipped
         */                      
-        while((SearchListAct!=NULL))// && counter<363)
+        while( (SearchListAct!=NULL))//&& counter<490)
         {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
                 fAct = SearchListAct;
@@ -3271 +3422 @@
 //                      gcTmp->getCodim2Normals(*gcTmp);
                         gcTmp->getExtremalRays(*gcTmp);
-                        gcTmp->normalize();// will be done in g2c 
+//                      gcTmp->normalize();// will be done in g2c 
                         //gcTmp->ddFacets should not be needed anymore, so
 //                      //NOTE If flip2 is used we need to get an interior point of gcTmp
@@ -3283 +3434 @@
 #endif
                         /*add facets to SLA here*/
+#ifdef SHALLOW
+                        SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
+#else
                         SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
+#endif
                         if(gfanHeuristic==1)
                         {
@@ -3377 +3532 @@
 /** \brief Make a set of rational vectors into integers
  *
- * We compute the lcm of the denominators and multiply with this to get integer values.         
+ * We compute the lcm of the denominators and multiply with this to get integer values.
+ * If the gcd of the nominators > 1 we divide by the gcd => primitive vector
  * \param dd_MatrixPtr,intvec
  */
@@ -3423 +3579 @@
                         
 //                      mpq_canonicalize(qkgV);
+        int ggT=1;
         for (int ii=0;ii<(M->colsize)-1;ii++)
         {
                 mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
                 n[ii]=(int)mpz_get_d(mpq_numref(res));
+                ggT=intgcd(ggT,n[ii]);
+        }
+        //Normalization
+        if(ggT>1)
+        {
+                for(int ii=0;ii<this->numVars;ii++)
+                        n[ii] /= ggT;
         }
         delete [] denom;
@@ -3541 +3705 @@
          */
         volatile bool removalOccured=FALSE;
-        int ctr=0;      //encountered equalities in SLA
-        int notParallelCtr=0;
+//      int ctr=0;      //encountered equalities in SLA
+//      int notParallelCtr=0;
 //      gcone::lengthOfSearchList=1;
         while(slEnd->next!=NULL)
@@ -3555 +3719 @@
                 {
                         intvec *fNormal=NULL;
-                        fNormal=fAct->getFacetNormal();
+                        fNormal=fAct->getFacetNormal();                 
                         slAct = slHead;
-                        notParallelCtr=0;
+                        //notParallelCtr=0;
                         /*If slAct==NULL and fAct!=NULL 
                         we just copy all remaining facets into SLA*/
@@ -3570 +3734 @@
                                                 if(slAct==NULL)
                                                 {
-                                                        slAct = new facet(*fCopy);//copy constructor
+                                                        slAct = new facet(*fCopy/*,TRUE*/);//copy constructor
                                                         slHead = slAct;                                                         
                                                 }
                                                 else
                                                 {
-                                                        slAct->next = new facet(*fCopy);
+                                                        slAct->next = new facet(*fCopy/*,TRUE*/);
                                                         slAct = slAct->next;
                                                 }
@@ -3592 +3756 @@
                                 cout << "Checking facet (";fNormal->show(1,1);cout << ") against (";slNormal->show(1,1);cout << ")" << endl;
 #endif                          
-                                if(areEqual(fAct,slAct))
-                                {
+//                              if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
+//                                      exit(-1);
+                                if(areEqual2(fAct,slAct))
+                                {                                       
                                         deleteMarker = slAct;
                                         if(slAct==slHead)
@@ -3648 +3814 @@
                                 facet *f2Act;
                                 f2Act = fAct->codim2Ptr;
-                                                
+
                                 slEnd->next = new facet();
                                 slEnd = slEnd->next;//Update slEnd
@@ -3663 +3829 @@
                                 //Better: have slEnd->interiorPoint point to fAct->interiorPoint
                                 //NOTE Only reference -> c.f. copy constructor
-//                              slEnd->setInteriorPoint(fAct->getInteriorPoint());
+                                //slEnd->setInteriorPoint(fAct->getInteriorPoint());
                                 slEnd->interiorPoint = fAct->interiorPoint;
                                 slEnd->prev = marker;
                                 //Copy codim2-facets
-//                              intvec *f2Normal=new intvec(this->numVars);
+                                //intvec *f2Normal=new intvec(this->numVars);
                                 while(f2Act!=NULL)
                                 {
@@ -3707 +3873 @@
         return slHead;
 }//addC2N
+
+/** Try using shallow copies*/
+facet * gcone::enqueue2(facet *f)
+{
+        assert(f!=NULL);
+#ifdef gfanp
+        timeval start, end;
+        gettimeofday(&start, 0);
+#endif
+        facet *slHead;
+        slHead = f;
+        facet *slAct;   //called with f=SearchListRoot
+        slAct = f;
+        facet *slEnd;   //Pointer to end of SLA
+        slEnd = f;
+        
+        facet *fAct;
+        fAct = this->facetPtr;//New facets to compare
+        facet *codim2Act;
+        codim2Act = this->facetPtr->codim2Ptr;
+//      facet *sl2Act;
+//      sl2Act = f->codim2Ptr;
+        volatile bool removalOccured=FALSE;
+        while(slEnd->next!=NULL)
+        {
+                slEnd=slEnd->next;
+        }       
+        while(fAct!=NULL)
+        {
+                if(fAct->isFlippable)
+                {
+                        slAct = slHead;
+                        if(slAct==NULL)
+                        {
+                                facet *fCopy;
+                                fCopy = fAct;                           
+                                while(fCopy!=NULL)
+                                {
+                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
+                                        {
+                                                if(slAct==NULL)
+                                                {
+                                                        slAct = slAct->shallowCopy(*fCopy);//shallow copy constructor
+                                                        slHead = slAct;                                                         
+                                                }
+                                                else
+                                                {
+                                                        slAct->next = slAct->shallowCopy(*fCopy);
+                                                        slAct = slAct->next;
+                                                }
+                                        }
+                                        fCopy = fCopy->next;
+                                }
+                                break;  
+                        }
+                        /*Comparison starts here*/
+                        while(slAct!=NULL)
+                        {
+                                removalOccured=FALSE;
+#ifdef gfan_DEBUG
+cout << "Checking facet (";fAct->fNormal->show(1,1);cout << ") against (";slAct->fNormal->show(1,1);cout << ")" << endl;
+#endif  
+                                if(areEqual2(fAct,slAct))
+                                {
+                                        facet *marker=slAct;
+                                        if(slAct==slHead)
+                                        {                                               
+                                                slHead = slAct->next;                                           
+                                                if(slHead!=NULL)
+                                                        slHead->prev = NULL;
+                                        }
+                                        else if (slAct==slEnd)
+                                        {
+                                                slEnd=slEnd->prev;
+                                                slEnd->next = NULL;
+                                        }                                                               
+                                        else
+                                        {
+                                                slAct->prev->next = slAct->next;
+                                                slAct->next->prev = slAct->prev;
+                                        }
+                                        removalOccured=TRUE;
+                                        gcone::lengthOfSearchList--;
+#ifdef gfan_DEBUG
+                                        cout << "Removing (";fAct->fNormal->show(1,1);cout << ") from list" << endl;
+#endif
+                                        marker->shallowDelete();
+//                                      delete(marker);
+                                        break;
+                                }
+                                slAct = slAct->next;
+                        }//while(slAct!=NULL)
+                        if(removalOccured==FALSE)
+                        {
+                                facet *marker=slEnd;
+                                slEnd->next = fAct->shallowCopy(*fAct);
+                                slEnd = slEnd->next;
+                                slEnd->prev=marker;
+                                gcone::lengthOfSearchList++;
+                        }
+                        fAct = fAct->next;
+                }
+                else
+                        fAct = fAct->next;
+        }//while(fAct!=NULL)
+        
+#ifdef gfanp
+        gettimeofday(&end, 0);
+        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
+#endif  
+        return slHead;
+}
+
 /** 
 * During flip2 every gc->baseRing gets two ringorder_a
@@ -3758 +4037 @@
         rChangeCurrRing(srcRing);                       
 }
-//A try with the STL
-facet * gcone::enqueue2(facet *f)
-{
-//      vector<facet*> searchList;
-//      facet *fAct;
-//      fAct=this->facetPtr;
-//      facet *slAct;
-//      slAct = f;
-//      facet *slHead;
-//      while(fAct!=NULL)
-//      {
-//              if(fAct->isFlippable==TRUE)
-//              {
-//                      if(searchList.size()==0)
-//                      {
-//                              facet *fCopy;
-//                              fCopy=fAct;
-//                              while(fCopy!=NULL)
-//                              {
-//                                      if(fCopy->isFlippable==TRUE)
-//                                              searchList.push_back(fCopy);
-//                                      fCopy=fCopy->next;
-//                              }
-//                      }
-//              }
-//      }
-}
 
 /** \brief Compute the gcd of two ints
  */
-inline int gcone::intgcd(const int a, const int b)
+inline int gcone::intgcd(const int &a, const int &b)
 {
         int r, p=a, q=b;
@@ -3869 +4121 @@
                                                 
         ofstream gcOutputFile(filename.c_str());
+        assert(gcOutputFile);
         facet *fAct;
         fAct = gc.facetPtr;                     
@@ -4167 +4420 @@
 }
         
-ring gcone::getBaseRing()
-{
-        return rCopy(this->baseRing);
-}
+
 /** \brief Sort the rays of a facet lexicographically
 */
@@ -4191 +4441 @@
 *
 */
-lists lprepareResult(gcone *gc, int n)
+lists lprepareResult(gcone *gc, const int n)
 {
         gcone *gcAct;
@@ -4205 +4455 @@
                 {
                         gcAct->readConeFromFile(gcAct->getUCN(),gcAct); 
-                        rChangeCurrRing(gcAct->getBaseRing());
+                        rChangeCurrRing(gcAct->getBaseRing());//NOTE memleak?
                 }               
                 res->m[ii].rtyp=LIST_CMD;
@@ -4224 +4474 @@
 //              rChangeCurrRing(tmpRing);
 //              l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
-                l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());
+                l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());//NOTE memleak?
 //              rChangeCurrRing(saveRing);
 
                 l->m[2].rtyp=INTVEC_CMD;
-                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
+                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));//NOTE memleak?
                 l->m[2].data=(void*)ivCopy(&iv);
                 
@@ -4311 +4561 @@
 float gcone::time_getCodim2Normals;
 float gcone::t_getExtremalRays;
+float gcone::t_ddPolyh;
 float gcone::time_flip;
 float gcone::time_flip2;
@@ -4322 +4573 @@
 float gcone::t_mI;
 float gcone::t_iP;
+float gcone::t_isParallel;
 unsigned gcone::parallelButNotEqual=0;
 unsigned gcone::numberOfFacetChecks=0;
@@ -4421 +4673 @@
 //      cout << "  t_iP:" << gcone::t_iP << endl;
         cout << "t_getExtremalRays:" << gcone::t_getExtremalRays << endl;
+        cout << "  t_ddPolyh:" << gcone::t_ddPolyh << endl;
         cout << "t_Flip:" << gcone::time_flip << endl;
         cout << "  t_markings:" << gcone::t_markings << endl;
@@ -4431 +4684 @@
         cout << "t_enqueue:" << gcone::time_enqueue << endl;
         cout << "  t_areEqual:" <<gcone::t_areEqual << endl;
+        cout << "t_isParallel:" <<gcone::t_isParallel << endl;
         cout << endl;
         cout << "Checked " << gcone::numberOfFacetChecks << " Facets" << endl;

kernel/gfan.h

@@ -23 +23 @@
 extern int gfanHeuristic;
 // extern dd_MatrixPtr ddLinealitySpace;
-#define gfanp
+
 // #ifdef gfanp
 // extern       static float time_getConeNormals;
@@ -75 +75 @@
                 facet *codim2Ptr;       //Pointer to (codim-2)-facet. Bit of recursion here ;-)
                 int numCodim2Facets;    //#of (codim-2)-facets of this facet. Set in getCodim2Normals()
+                unsigned numRays;       //Number of spanning rays of the cone
                 ring flipRing;          //the ring on the other side of the facet
-                unsigned numRays;       //Number of spanning rays of the cone                   
+                                        
                 /** The default constructor. */
                 facet();
                 /** Constructor for lower dimensional faces*/
-                facet(int const &n);
+                facet(const int &n);
                 /**  The copy constructor */
                 facet(const facet& f);
-                facet(const facet& f, bool shallow);
+                facet* shallowCopy(const facet& f);
+                void shallowDelete();
                 /** The default destructor */
                 ~facet();
@@ -113 +115 @@
                 inline void setInteriorPoint(intvec *iv);
                 inline intvec *getInteriorPoint();
+                inline const intvec *getRef2InteriorPoint();
                 /** \brief Debugging function
                  * prints the facet normal an all (codim-2)-facets that belong to it
@@ -145 +148 @@
                 static float time_getCodim2Normals;
                 static float t_getExtremalRays;
+                static float t_ddPolyh;
                 static float time_flip;
                 static float time_flip2;
@@ -157 +161 @@
                 static float t_mI;
                 static float t_iP;
+                static float t_isParallel;
                 static unsigned parallelButNotEqual;
                 static unsigned numberOfFacetChecks;
@@ -196 +201 @@
                 inline ring getBaseRing();
                 inline void setIntPoint(intvec *iv);
-                inline intvec *getIntPoint();
+                inline intvec *getIntPoint(bool shallow=FALSE);
                 inline void showIntPoint();
                 inline void setNumFacets();
@@ -213 +218 @@
 //              inline int dotProduct(const intvec* a, const intvec *b);
 //              inline bool isParallel(const intvec* a, const intvec* b);
-                inline bool areEqual(intvec &a, intvec &b);
+                inline bool ivAreEqual(intvec &a, intvec &b);
                 inline bool areEqual( facet *f,  facet *g);
                 inline bool areEqual2(facet* f, facet *g);
-                inline int intgcd(int a, int b);
+                inline int intgcd(const int &a, const int &b);
                 inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE);
                 inline void readConeFromFile(int gcNum, gcone *gc);
@@ -252 +257 @@
                 friend class facet;     
 };
-lists lprepareResult(gcone *gc, int n);
+lists lprepareResult(gcone *gc, const int n);
 // bool iv64isStrictlyPositive(intvec *);
 #endif