Index: kernel/gfan.cc
===================================================================
--- kernel/gfan.cc	(revision c57a83edd49526195f04ddf602fca071cafed3b7)
+++ kernel/gfan.cc	(revision 73809b3238e9c892b9499243f1d22f8a468631a8)
@@ -75,7 +75,13 @@
 #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,5 +120,5 @@
 * 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,5 +150,5 @@
 	//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,5 +157,5 @@
 	while(f2Copy!=NULL)
 	{
-		if(f2Act==NULL)
+		if(f2Act==NULL || f2Act==(facet*)0xfefefefefefefefe)
 		{
 			f2Act=new facet(2);
@@ -175,7 +181,31 @@
 
 /** \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,6 +247,7 @@
 	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,27 +256,41 @@
 	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,5 +430,5 @@
 // 	return this->fNormal;
 }
-		
+
 /** Method to print the facet normal*/
 inline void facet::printNormal()
@@ -456,5 +501,10 @@
 {
 	return ivCopy(this->interiorPoint);
-}	
+}
+
+inline const intvec *facet::getRef2InteriorPoint()
+{
+	return (this->interiorPoint);
+}
 		
 /** \brief Debugging function
@@ -565,6 +615,17 @@
 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,5 +650,6 @@
 	//dd_FreeMatrix(this->ddFacets);
 }			
-	
+
+/** Returns the number of cones existing at the time*/
 inline int gcone::getCounter()
 {
@@ -603,8 +665,11 @@
 }
 		
-/** \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,4 +827,9 @@
 {
 	return this->pred;
+}
+/** Returns a copy of the this->baseRing */
+inline ring gcone::getBaseRing()
+{
+	return rCopy(this->baseRing);
 }
 
@@ -796,6 +866,7 @@
 	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,4 +1076,5 @@
 	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,7 +1082,12 @@
 			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,4 +1389,5 @@
 * 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,4 +1396,6 @@
 	timeval start, end;
 	gettimeofday(&start, 0);
+	timeval poly_start, poly_end;
+	gettimeofday(&poly_start,0);
 #endif
 	//Add lineality space - dd_LinealitySpace
@@ -1347,11 +1427,17 @@
 	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,4 +1458,5 @@
 	}
 	//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,4 +1476,5 @@
 	mpz_clear(den);
 	mpz_clear(tmp);
+	int ggT=(*ivIntPointOfCone)[0];
 	for (int ii=0;ii<(this->numVars);ii++)
 	{
@@ -1396,6 +1484,15 @@
 		(*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,6 +1510,15 @@
 		}
 		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,5 +1536,5 @@
 		{			
 			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,5 +1562,6 @@
 			}
 			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,9 +1584,17 @@
 				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,6 +1752,6 @@
 	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,5 +2460,5 @@
 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,5 +2475,9 @@
  */
 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,4 +2493,8 @@
 		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,4 +3075,5 @@
 }
 		
+//NOTE not needed anywhere
 ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
 {
@@ -2973,5 +3097,5 @@
 /** \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,8 +3226,4 @@
 }	
 
-/** 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,5 +3274,5 @@
 				
 	//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,4 +3287,21 @@
 		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,5 +3354,6 @@
 			}
 			fAct = fAct->next;
-			delete fNormal;			
+			delete fNormal;
+#endif
 		}//if(fAct->isFlippable==TRUE)			
 		else {fAct = fAct->next;}
@@ -3224,8 +3362,21 @@
 	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,5 +3394,5 @@
 	  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,5 +3422,5 @@
 // 			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,5 +3434,9 @@
 #endif
 			/*add facets to SLA here*/
+#ifdef SHALLOW
+			SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
+#else
 			SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
+#endif
  			if(gfanHeuristic==1)
  			{
@@ -3377,5 +3532,6 @@
 /** \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,8 +3579,16 @@
 			
 // 			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,6 +3705,6 @@
 	 */
 	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,7 +3719,7 @@
 		{
 			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,10 +3734,10 @@
 						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,6 +3756,8 @@
 				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,5 +3814,5 @@
 				facet *f2Act;
 				f2Act = fAct->codim2Ptr;
-						
+
 				slEnd->next = new facet();
 				slEnd = slEnd->next;//Update slEnd
@@ -3663,9 +3829,9 @@
 				//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,4 +3873,117 @@
 	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,35 +4037,8 @@
 	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,4 +4121,5 @@
 						
 	ofstream gcOutputFile(filename.c_str());
+	assert(gcOutputFile);
 	facet *fAct;
 	fAct = gc.facetPtr;			
@@ -4167,8 +4420,5 @@
 }
 	
-ring gcone::getBaseRing()
-{
-	return rCopy(this->baseRing);
-}
+
 /** \brief Sort the rays of a facet lexicographically
 */
@@ -4191,5 +4441,5 @@
 *
 */
-lists lprepareResult(gcone *gc, int n)
+lists lprepareResult(gcone *gc, const int n)
 {
 	gcone *gcAct;
@@ -4205,5 +4455,5 @@
 		{
 			gcAct->readConeFromFile(gcAct->getUCN(),gcAct);	
-			rChangeCurrRing(gcAct->getBaseRing());
+			rChangeCurrRing(gcAct->getBaseRing());//NOTE memleak?
 		}		
 		res->m[ii].rtyp=LIST_CMD;
@@ -4224,9 +4474,9 @@
 // 		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,4 +4561,5 @@
 float gcone::time_getCodim2Normals;
 float gcone::t_getExtremalRays;
+float gcone::t_ddPolyh;
 float gcone::time_flip;
 float gcone::time_flip2;
@@ -4322,4 +4573,5 @@
 float gcone::t_mI;
 float gcone::t_iP;
+float gcone::t_isParallel;
 unsigned gcone::parallelButNotEqual=0;
 unsigned gcone::numberOfFacetChecks=0;
@@ -4421,4 +4673,5 @@
 // 	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,4 +4684,5 @@
 	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;
Index: kernel/gfan.h
===================================================================
--- kernel/gfan.h	(revision c57a83edd49526195f04ddf602fca071cafed3b7)
+++ kernel/gfan.h	(revision 73809b3238e9c892b9499243f1d22f8a468631a8)
@@ -23,5 +23,5 @@
 extern int gfanHeuristic;
 // extern dd_MatrixPtr ddLinealitySpace;
-#define gfanp
+
 // #ifdef gfanp
 // extern	static float time_getConeNormals;
@@ -75,13 +75,15 @@
 		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,4 +115,5 @@
 		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,4 +148,5 @@
 		static float time_getCodim2Normals;
 		static float t_getExtremalRays;
+		static float t_ddPolyh;
 		static float time_flip;
 		static float time_flip2;
@@ -157,4 +161,5 @@
 		static float t_mI;
 		static float t_iP;
+		static float t_isParallel;
 		static unsigned parallelButNotEqual;
 		static unsigned numberOfFacetChecks;
@@ -196,5 +201,5 @@
 		inline ring getBaseRing();
 		inline void setIntPoint(intvec *iv);
-		inline intvec *getIntPoint();
+		inline intvec *getIntPoint(bool shallow=FALSE);
 		inline void showIntPoint();
 		inline void setNumFacets();
@@ -213,8 +218,8 @@
 // 		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,5 +257,5 @@
 		friend class facet;	
 };
-lists lprepareResult(gcone *gc, int n);
+lists lprepareResult(gcone *gc, const int n);
 // bool iv64isStrictlyPositive(intvec *);
 #endif
