Changeset 4199b3 in git

Timestamp:

Mar 17, 2011, 5:37:55 PM (12 years ago)

Author:

Martin Monerjan

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')

Children:

e7463f6a2c8fb0cc7abe50aa54f387c997fec2d8

Parents:

9a038bd3951d010760b79ca2e8b0eb43932c0b9a

Message:

HAVE_GFAN merged with HAVE_FANS
defined USE_ZFAN
gcRays2IntMat
prepareGfanLib
cleanup


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

Location:

kernel

Files:

• gfan.cc (modified) (236 diffs)
• gfan.h (modified) (16 diffs)

kernel/gfan.cc

@@ -9 +9 @@
 #include <kernel/mod2.h>
 
-#ifdef HAVE_GFAN
+#ifdef HAVE_FANS
 #include <kernel/options.h>
 #include <kernel/kstd1.h>
 #include <kernel/kutil.h>
- #include "int64vec.h"
+// #include "int64vec.h"
 #include <kernel/polys.h>
 #include <kernel/ideals.h>
 #include <kernel/kmatrix.h>
 #include <kernel/GMPrat.h>
-#include "fast_maps.h"  //Mapping of ideals
- #include "maps.h"
-#include "ring.h"       apparently not needed
- #include "structs.h"
+
+#include "ring.h"       //apparently not needed
+// #include "structs.h"
 #include <Singular/lists.h>
 #include <kernel/prCopy.h>
 #include <kernel/stairc.h>
- #include <bitset>
-#include <fstream>      read-write cones to files
- #include <gmp.h>
+// #include <bitset>
+#include <fstream>      //read-write cones to files
+// #include <gmp.h>
 #include <string>
 #include <sstream>
- #include <time.h>
+// #include <time.h>
 #include <stdlib.h>
-#include <gmpxx.h>
- #include <vector>
 #include <assert.h>
-
+// #include <Singular/bbfan.h>
+// #include <Singular/bbcone.h>
+#include <gfanlib/gfanlib.h>
 
 /*DO NOT REMOVE THIS*/
@@ -42 +41 @@
 #endif
 
-Hacks for different working places
+//Hacks for different working places
 #define p800
-
-#ifdef UNI
-#include <ers/urmel/alggeom/monerjan/cddlib/include/setoper.h>
-#include <ers/urmel/alggeom/monerjan/cddlib/include/cdd.h>
-#endif
-
-#ifdef HOME
-#include <me/momo/studium/diplomarbeit/cddlib/include/setoper.h>
-#include <me/momo/studium/diplomarbeit/cddlib/include/cdd.h>
-#endif
-
-#ifdef ITWM
-#include <slg/monerjan/cddlib/include/setoper.h>
-#include <slg/monerjan/cddlib/include/cdd.h>
-#include <slg/monerjan/cddlib/include/cddmp.h>
-#endif
 
 #ifdef p800
@@ -68 +51 @@
 
 #ifndef gfan_DEBUG
- #define gfan_DEBUG
+// #define gfan_DEBUG
 #ifndef gfan_DEBUGLEVEL
 #define gfan_DEBUGLEVEL 1
@@ -74 +57 @@
 #endif
 
-NOTE Defining this will slow things down!
-Only good for very coarse profiling
- #define gfanp
+//NOTE Defining this will slow things down!
+//Only good for very coarse profiling
+// #define gfanp
 #ifdef gfanp
 #include <sys/time.h>
@@ -84 +67 @@
 #ifndef SHALLOW
 #define SHALLOW
+#endif
+
+#ifndef USE_ZFAN
+#define USE_ZFAN
 #endif
 
@@ -96 +83 @@
 *
 */
-
+                
 /** \brief The default constructor for facets
 */
-facet::facet()
-{
-         Pointer to next facet.  */
+facet::facet()                  
+{
+        // Pointer to next facet.  */
         /* Defaults to NULL. This way there is no need to check explicitly */
         this->fNormal=NULL;
-        this->interiorPoint=NULL;
+        this->interiorPoint=NULL;               
         this->UCN=0;
         this->codim2Ptr=NULL;
-        this->codim=1;          default for (codim-1)-facets
+        this->codim=1;          //default for (codim-1)-facets
         this->numCodim2Facets=0;
         this->numRays=0;
         this->flipGB=NULL;
         this->next=NULL;
-        this->prev=NULL;
-        this->flipRing=NULL;    the ring on the other side
+        this->prev=NULL;                
+        this->flipRing=NULL;    //the ring on the other side
         this->isFlippable=FALSE;
 }
-
+                
 /** \brief Constructor for facets of codim == 2
 * Note that as of now the code of the constructors is only for facets and codim2-faces. One
@@ -124 +111 @@
 {
         this->fNormal=NULL;
-        this->interiorPoint=NULL;
+        this->interiorPoint=NULL;                       
         this->UCN=0;
         this->codim2Ptr=NULL;
@@ -130 +117 @@
         {
                 this->codim=n;
-        }NOTE Handle exception here!
+        }//NOTE Handle exception here!                  
         this->numCodim2Facets=0;
         this->numRays=0;
@@ -139 +126 @@
         this->isFlippable=FALSE;
 }
-
+                
 /** \brief The copy constructor
 * By default only copies the fNormal, f2Normals and UCN
@@ -148 +135 @@
         this->UCN=f.UCN;
         this->isFlippable=f.isFlippable;
-        Needed for flip2
-        NOTE ONLY REFERENCE
-        this->interiorPoint=iv64Copy(f.interiorPoint);only referencing is prolly not the best thing to do in a copy constructor
+        //Needed for flip2
+        //NOTE ONLY REFERENCE
+        this->interiorPoint=iv64Copy(f.interiorPoint);//only referencing is prolly not the best thing to do in a copy constructor
         facet* f2Copy;
         f2Copy=f.codim2Ptr;
@@ -168 +155 @@
                 {
                         f2Act=new facet(2);
-                        this->codim2Ptr=f2Act;
+                        this->codim2Ptr=f2Act;                  
                 }
                 else
@@ -180 +167 @@
                 int64vec *f2Normal;
                 f2Normal = f2Copy->getFacetNormal();
-                f2Act->setFacetNormal(f2Copy->getFacetNormal());
+//              f2Act->setFacetNormal(f2Copy->getFacetNormal());
                 f2Act->setFacetNormal(f2Normal);
                 delete f2Normal;
                 f2Act->setUCN(f2Copy->getUCN());
                 f2Copy = f2Copy->next;
-        }
+        }       
 }
 
@@ -209 +196 @@
 {
 #ifdef gfan_DEBUG
-        printf("shallowdel@UCN %i\n", this->getUCN());
+//      printf("shallowdel@UCN %i\n", this->getUCN());
 #endif
         this->fNormal=NULL;
-        this->UCN=0;
+//      this->UCN=0;
         this->interiorPoint=NULL;
         this->codim2Ptr=NULL;
@@ -219 +206 @@
         this->flipGB=NULL;
         this->flipRing=NULL;
-        assert(this->fNormal==NULL);
-        delete(this);
-}
-
+        assert(this->fNormal==NULL);    
+//      delete(this);
+}
+                
 /** The default destructor */
 facet::~facet()
 {
 #ifdef gfan_DEBUG
-        printf("~facet@UCN %i\n",this->getUCN());
+//      printf("~facet@UCN %i\n",this->getUCN());
 #endif
         if(this->fNormal!=NULL)
@@ -234 +221 @@
                 delete this->interiorPoint;
         /* Cleanup the codim2-structure */
-        if(this->codim==2)
-        {
-                facet *codim2Ptr;
-                codim2Ptr = this->codim2Ptr;
-                while(codim2Ptr!=NULL)
-                {
-                        if(codim2Ptr->fNormal!=NULL)
-                        {
-                                delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
-                                codim2Ptr = codim2Ptr->next;
-                        }
-                }
-        }
-        The rays are stored in the cone!
+//      if(this->codim==2)
+//      {
+//              facet *codim2Ptr;
+//              codim2Ptr = this->codim2Ptr;
+//              while(codim2Ptr!=NULL)
+//              {
+//                      if(codim2Ptr->fNormal!=NULL)
+//                      {
+//                              delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
+//                              codim2Ptr = codim2Ptr->next;
+//                      }
+//              }
+//      }
+        //The rays are stored in the cone!
         if(this->flipGB!=NULL)
                 idDelete((ideal *)&this->flipGB);
-        if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
-                rDelete(this->flipRing); //See vol II/134
-        this->flipRing=NULL;
+//      if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
+//              rDelete(this->flipRing); //See vol II/134
+//      this->flipRing=NULL;
         this->prev=NULL;
         this->next=NULL;
 }
-
+                
 inline const int64vec *facet::getRef2FacetNormal() const
 {
         return(this->fNormal);
-}
+}       
 
 /** Equality check for facets based on unique interior points*/
@@ -281 +268 @@
                 }
         }
-        if( fIntP->compare(gIntP)!=0) res=FALSE;
+//      if( fIntP->compare(gIntP)!=0) res=FALSE;
 #ifdef gfanp
         gettimeofday(&end, 0);
         gcone::t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
-#endif
+#endif  
         return res;
 }
@@ -308 +295 @@
         int notParallelCtr=0;
         int ctr=0;
-        const int64vec* fNormal; No new since iv64Copy and therefore getFacetNormal return a new
+        const int64vec* fNormal; //No new since iv64Copy and therefore getFacetNormal return a new
         const int64vec* sNormal;
-        fNormal = f->getRef2FacetNormal();->getFacetNormal();
-        sNormal = s->getRef2FacetNormal();->getFacetNormal();
-        Do not need parallelity. Too time consuming
-        if(!isParallel(*fNormal,*sNormal))
-        if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
-                notParallelCtr++;
-        else//parallelity, so we check the codim2-facets
+        fNormal = f->getRef2FacetNormal();//->getFacetNormal();
+        sNormal = s->getRef2FacetNormal();//->getFacetNormal();
+#include "intvec.h"
+        //Do not need parallelity. Too time consuming
+//      if(!isParallel(*fNormal,*sNormal))
+//      if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
+ //             notParallelCtr++;
+//      else//parallelity, so we check the codim2-facets
         int64vec *fNRef=const_cast<int64vec*>(fNormal);
         int64vec *sNRef=const_cast<int64vec*>(sNormal);
-        if(isParallel(*fNormal,*sNormal))
+//      if(isParallel(*fNormal,*sNormal))       
         if(isParallel(*fNRef,*sNRef))
-        if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
+//      if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
         {
                 facet* f2Act;
                 facet* s2Act;
-                f2Act = f->codim2Ptr;
+                f2Act = f->codim2Ptr;           
                 ctr=0;
                 while(f2Act!=NULL)
                 {
                         const int64vec* f2Normal;
-                        f2Normal = f2Act->getRef2FacetNormal();->getFacetNormal();
-                        int64vec *f2Ref=const_cast<int64vec*>(f2Normal);
+                        f2Normal = f2Act->getRef2FacetNormal();//->getFacetNormal();
+//                      int64vec *f2Ref=const_cast<int64vec*>(f2Normal);
                         s2Act = s->codim2Ptr;
                         while(s2Act!=NULL)
                         {
                                 const int64vec* s2Normal;
-                                s2Normal = s2Act->getRef2FacetNormal();->getFacetNormal();
-                                bool foo=areEqual(f2Normal,s2Normal);
-                                int64vec *s2Ref=const_cast<int64vec*>(s2Normal);
+                                s2Normal = s2Act->getRef2FacetNormal();//->getFacetNormal();
+//                              bool foo=areEqual(f2Normal,s2Normal);
+//                              int64vec *s2Ref=const_cast<int64vec*>(s2Normal);
                                 int foo=f2Normal->compare(s2Normal);
                                 if(foo==0)
                                         ctr++;
                                 s2Act = s2Act->next;
-                                delete s2Normal;
+//                              delete s2Normal;
                         }
-                        delete f2Normal;
+//                      delete f2Normal;
                         f2Act = f2Act->next;
-                }
-        }
-        delete fNormal;
-        delete sNormal;
+                }               
+        }
+//      delete fNormal;
+//      delete sNormal;
         if(ctr==f->numCodim2Facets)
                 res=TRUE;
@@ -365 +353 @@
 #endif
         return res;
-        int64vec *foo=ivNeg(sNormal);
-        if(fNormal->compare(foo)!=0) //facet normals
-        {
-                delete foo;
-                res=FALSE;
-        }
-        else
-        {
-                facet* f2Act;
-                facet* s2Act;
-                f2Act = f->codim2Ptr;
-                ctr=0;
-                while(f2Act!=NULL)
-                {
-                        int64vec* f2Normal;
-                        f2Normal = f2Act->getFacetNormal();
-                        s2Act = s->codim2Ptr;
-                        while(s2Act!=NULL)
-                        {
-                                int64vec* s2Normal;
-                                s2Normal = s2Act->getFacetNormal();
- //                             bool foo=areEqual(f2Normal,s2Normal);
-                                int foo=f2Normal->compare(s2Normal);
-                                if(foo==0)
-                                        ctr++;
-                                s2Act = s2Act->next;
-                                delete s2Normal;
-                        }
-                        delete f2Normal;
-                        f2Act = f2Act->next;
-                }
-        }
-        delete fNormal;
-        delete sNormal;
-        if(ctr==f->numCodim2Facets)
-                res=TRUE;
-        return res;
-}
-
+//      int64vec *foo=ivNeg(sNormal);
+//      if(fNormal->compare(foo)!=0) //facet normals
+//      {
+//              delete foo;
+//              res=FALSE;
+//      }
+//      else
+//      {
+//              facet* f2Act;
+//              facet* s2Act;
+//              f2Act = f->codim2Ptr;           
+//              ctr=0;
+//              while(f2Act!=NULL)
+//              {
+//                      int64vec* f2Normal;
+//                      f2Normal = f2Act->getFacetNormal();
+//                      s2Act = s->codim2Ptr;
+//                      while(s2Act!=NULL)
+//                      {
+//                              int64vec* s2Normal;
+//                              s2Normal = s2Act->getFacetNormal();
+// //                           bool foo=areEqual(f2Normal,s2Normal);
+//                              int foo=f2Normal->compare(s2Normal);
+//                              if(foo==0)
+//                                      ctr++;
+//                              s2Act = s2Act->next;
+//                              delete s2Normal;
+//                      }
+//                      delete f2Normal;
+//                      f2Act = f2Act->next;
+//              }
+//      }               
+//      delete fNormal;
+//      delete sNormal;
+//      if(ctr==f->numCodim2Facets)
+//              res=TRUE;
+//      return res;
+}       
+                
 /** Stores the facet normal \param int64vec*/
 inline void facet::setFacetNormal(int64vec *iv)
@@ -409 +397 @@
         if(this->fNormal!=NULL)
                 delete this->fNormal;
-        this->fNormal = iv64Copy(iv);
-}
-
-/** Hopefully returns the facet normal
+        this->fNormal = iv64Copy(iv);                   
+}
+                
+/** Hopefully returns the facet normal 
 * Mind: iv64Copy returns a new int64vec, so use this in the following way:
 * int64vec *iv;
@@ -420 +408 @@
 */
 inline int64vec *facet::getFacetNormal() const
-{
+{                               
         return iv64Copy(this->fNormal);
-        return this->fNormal;
+//      return this->fNormal;
 }
 
@@ -430 +418 @@
         fNormal->show();
 }
-
+                
 /** Store the flipped GB*/
 inline void facet::setFlipGB(ideal I)
@@ -436 +424 @@
         this->flipGB=idCopy(I);
 }
-
+                
 /** Returns a pointer to the flipped GB
 Seems not be used
@@ -445 +433 @@
         return this->flipGB;
 }
-
+                
 /** Print the flipped GB*/
 inline void facet::printFlipGB()
@@ -453 +441 @@
 #endif
 }
-
+                
 /** Set the UCN */
 inline void facet::setUCN(int n)
@@ -459 +447 @@
         this->UCN=n;
 }
-
-/** \brief Get the UCN
+                
+/** \brief Get the UCN 
 * Returns the UCN iff this != NULL, else -1
 */
@@ -474 +462 @@
 #ifdef NDEBUG
         if(this!=NULL)
-#endif
+#endif                  
                 return this->UCN;
         else
                 return -1;
 }
-
+                
 /** Store an interior point of the facet */
 inline void facet::setInteriorPoint(int64vec *iv)
@@ -487 +475 @@
         this->interiorPoint = iv64Copy(iv);
 }
-
+                
 /** Returns a pointer to this->interiorPoint
 * MIND: iv64Copy returns a new int64vec
@@ -501 +489 @@
         return (this->interiorPoint);
 }
-
+                
 /** \brief Debugging function
 * prints the facet normal an all (codim-2)-facets that belong to it
 */
 volatile void facet::fDebugPrint()
-{
-        facet *codim2Act;
+{                       
+        facet *codim2Act;                       
         codim2Act = this->codim2Ptr;
         int64vec *fNormal;
-        fNormal = this->getFacetNormal();
+        fNormal = this->getFacetNormal();       
         printf("=======================\n");
         printf("Facet normal = (");fNormal->show(1,1);printf(")\n");
@@ -524 +512 @@
         }
         printf("=======================\n");
-        delete fNormal;
-}
-
-friend class gcone;     //Bad style
+        delete fNormal; 
+}               
+                
+//friend class gcone;   //Bad style
 
 
@@ -538 +526 @@
 
 
-/** \brief Default constructor.
+/** \brief Default constructor. 
  *
  * Initialises this->next=NULL and this->facetPtr=NULL
@@ -546 +534 @@
         this->next=NULL;
         this->prev=NULL;
-        this->facetPtr=NULL;    maybe this->facetPtr = new facet();
+        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
         this->baseRing=currRing;
         this->counter++;
-        this->UCN=this->counter;
+        this->UCN=this->counter;                        
         this->numFacets=0;
         this->ivIntPt=NULL;
         this->gcRays=NULL;
 }
-
+                
 /** \brief Constructor with ring and ideal
  *
- * This constructor takes the root ring and the root ideal as parameters and stores
+ * This constructor takes the root ring and the root ideal as parameters and stores 
  * them in the private members gcone::rootRing and gcone::inputIdeal
  * This constructor is only called once in the computation of the Gröbner fan,
@@ -569 +557 @@
         this->next=NULL;
         this->prev=NULL;
-        this->facetPtr=NULL;
+        this->facetPtr=NULL;                    
         this->inputIdeal=I;
         this->baseRing=currRing;
@@ -579 +567 @@
         this->gcRays=NULL;
 }
-
-/** \brief Copy constructor
+                
+/** \brief Copy constructor 
  *
  * Copies a cone, sets this->gcBasis to the flipped GB
  * Call this only after a successful call to gcone::flip which sets facet::flipGB
-*/
+*/              
 gcone::gcone(const gcone& gc, const facet &f)
 {
         this->next=NULL;
-        this->prev=(gcone *)&gc; //comment in to get a tree
+//      this->prev=(gcone *)&gc; //comment in to get a tree
         this->prev=NULL;
-        this->numVars=gc.numVars;
+        this->numVars=gc.numVars;                                               
         this->counter++;
         this->UCN=this->counter;
@@ -596 +584 @@
         this->facetPtr=NULL;
         this->gcBasis=idrCopyR(f.flipGB, f.flipRing);
-        this->inputIdeal=idCopy(this->gcBasis);
+//      this->inputIdeal=idCopy(this->gcBasis);
         this->baseRing=rCopy(f.flipRing);
         this->numFacets=0;
@@ -602 +590 @@
         this->gcRays=NULL;
 }
-
-/** \brief Default destructor
+                
+/** \brief Default destructor 
 */
 gcone::~gcone()
@@ -619 +607 @@
                 idDelete((ideal *)&this->gcBasis);
 #endif
-        idDelete((ideal *)&this->gcBasis);
-        if(this->inputIdeal!=NULL)
-                idDelete((ideal *)&this->inputIdeal);
-        if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
-                rDelete(this->rootRing);
+//      idDelete((ideal *)&this->gcBasis);
+//      if(this->inputIdeal!=NULL)
+//              idDelete((ideal *)&this->inputIdeal);
+//      if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
+//              rDelete(this->rootRing);
         if(this->UCN!=1 && this->baseRing!=NULL)
                 rDelete(this->baseRing);
@@ -638 +626 @@
         }
         this->counter--;
-        should be deleted in noRevS
-        dd_FreeMatrix(this->ddFacets);
-        dd_FreeMatrix(this->ddFacets);
+        //should be deleted in noRevS
+//      dd_FreeMatrix(this->ddFacets);
+        //dd_FreeMatrix(this->ddFacets);
         for(int ii=0;ii<this->numRays;ii++)
                 delete(gcRays[ii]);
         omFree(gcRays);
-}
+}                       
 
 /** Returns the number of cones existing at the time*/
@@ -651 +639 @@
         return this->counter;
 }
-
+        
 /** \brief Set the interior point of a cone */
 inline void gcone::setIntPoint(int64vec *iv)
@@ -659 +647 @@
         this->ivIntPt=iv64Copy(iv);
 }
-
+                
 /** \brief Returns either a physical copy the interior point of a cone or just a reference to it.*/
 inline int64vec *gcone::getIntPoint(bool shallow)
@@ -668 +656 @@
                 return iv64Copy(this->ivIntPt);
 }
-
+                
 /** \brief Print the interior point */
 inline void gcone::showIntPoint()
@@ -674 +662 @@
         ivIntPt->show();
 }
-
+                
 /** \brief Print facets
  * This is mainly for debugging purposes. Usually called from within gdb
@@ -708 +696 @@
         printf("\n");
 }
-
+                
 /** For debugging purposes only */
 static volatile void showSLA(facet &f)
@@ -718 +706 @@
                 facet *codim2Act;
                 codim2Act = fAct->codim2Ptr;
-
+                
                 printf("\n");
                 while(fAct!=NULL)
                 {
-                        int64vec *fNormal;
+                        int64vec *fNormal;              
                         fNormal=fAct->getFacetNormal();
                         printf("(");fNormal->show(1,0);
@@ -745 +733 @@
         }
 }
-
+                
 static void idDebugPrint(const ideal &I)
 {
@@ -761 +749 @@
 static void invPrint(const ideal &I)
 {
-        int numElts=IDELEMS(I);
-        cout << "inv = ";
-        for(int ii=0;ii<numElts;ii++);
-        {
-                pWrite0(pHead(I->m[ii]));
-                cout << ",";
-        }
-        cout << endl;
+//      int numElts=IDELEMS(I);
+//      cout << "inv = ";
+//      for(int ii=0;ii<numElts;ii++);
+//      {
+//              pWrite0(pHead(I->m[ii]));
+//              cout << ",";
+//      }
+//      cout << endl;
 }
 
@@ -784 +772 @@
         return res;
 }
-
+                
 /** \brief Set gcone::numFacets */
 inline void gcone::setNumFacets()
 {
 }
-
+                
 /** \brief Get gcone::numFacets */
 inline int gcone::getNumFacets()
@@ -795 +783 @@
         return this->numFacets;
 }
-
+                
 inline int gcone::getUCN()
 {
-        if( this!=NULL) && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
+        if( this!=NULL)// && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
                 return this->UCN;
         else
@@ -829 +817 @@
  * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
  * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
- * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
+ * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects 
  * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
  * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
@@ -844 +832 @@
 #endif
         poly aktpoly;
-        int rows;                        will contain the dimensions of the ineq matrix - deprecated by
+        int rows;                       // will contain the dimensions of the ineq matrix - deprecated by
         dd_rowrange ddrows;
         dd_colrange ddcols;
-        dd_rowset ddredrows;             # of redundant rows in ddineq
-        dd_rowset ddlinset;              the opposite
-        dd_rowindex ddnewpos=NULL;                 all to make dd_Canonicalize happy
-        dd_NumberType ddnumb=dd_Integer; Number type
-        dd_ErrorType dderr=dd_NoError;
-        Compute the # inequalities i.e. rows of the matrix
-        rows=0; Initialization
+        dd_rowset ddredrows;            // # of redundant rows in ddineq
+        dd_rowset ddlinset;             // the opposite
+        dd_rowindex ddnewpos=NULL;                // all to make dd_Canonicalize happy
+        dd_NumberType ddnumb=dd_Integer; //Number type
+        dd_ErrorType dderr=dd_NoError;          
+        //Compute the # inequalities i.e. rows of the matrix
+        rows=0; //Initialization
         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;
         }
 
-        dd_rowrange aktmatrixrow=0;      needed to store the diffs of the expvects in the rows of ddineq
+        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
         ddrows=rows;
         ddcols=this->numVars;
-        dd_MatrixPtr ddineq;            Matrix to store the inequalities
-        ddineq=dd_CreateMatrix(ddrows,ddcols+1); The first col has to be 0 since cddlib checks for additive consts there
-
-         We loop through each g\in GB and compute the resulting inequalities
+        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                      
+        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
+                
+        // We loop through each g\in GB and compute the resulting inequalities
         for (int i=0; i<IDELEMS(I); i++)
         {
-                aktpoly=(poly)I->m[i];          get aktpoly as i-th component of I
-                simpler version of storing expvect diffs
+                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
+                //simpler version of storing expvect diffs
                 int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
                 pGetExpV(aktpoly,leadexpv);
@@ -879 +867 @@
                         pNextTerm/*aktpoly*/=pNext(pNextTerm);
                         int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
-                        pGetExpV(pNextTerm,tailexpv);
+                        pGetExpV(pNextTerm,tailexpv);                   
                         for(int kk=1;kk<=this->numVars;kk++)
-                        {
+                        {                               
                                 dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
                         }
                         aktmatrixrow += 1;
                         omFree(tailexpv);
-                }
-                omFree(leadexpv);
-        } for
+                }               
+                omFree(leadexpv);       
+        } //for
 #if true
         /*Let's make the preprocessing here. This could already be done in the above for-loop,
@@ -895 +883 @@
         * Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
         */
-        ideal initialForm=idInit(IDELEMS(I),1);
-        int64vec *gamma=new int64vec(this->numVars);
+//      ideal initialForm=idInit(IDELEMS(I),1);
+//      int64vec *gamma=new int64vec(this->numVars);
         int falseGammaCounter=0;
         int *redRowsArray=NULL;
         int num_alloc=0;
-        int num_elts=0;
+        int num_elts=0; 
         for(int ii=0;ii<ddineq->rowsize;ii++)
         {
                 ideal initialForm=idInit(IDELEMS(I),I->rank);
-                read row ii into gamma
-                int64 tmp;
+                //read row ii into gamma
+//              int64 tmp;
                 int64vec *gamma=new int64vec(this->numVars);
                 for(int jj=1;jj<=this->numVars;jj++)
@@ -915 +903 @@
                 computeInv((ideal&)I,initialForm,*gamma);
                 delete gamma;
-                Create leading ideal
+                //Create leading ideal
                 ideal L=idInit(IDELEMS(initialForm),1);
                 for(int jj=0;jj<IDELEMS(initialForm);jj++)
@@ -922 +910 @@
                         L->m[jj]=pCopy(/*pHead(initialForm->m[jj]))*/p);
                         pDelete(&p);
-                }
-
+                }               
+                
                 LObject *P = new sLObject();//TODO What's the difference between sLObject and LObject?
                 memset(P,0,sizeof(LObject));
@@ -930 +918 @@
                 {
                         bool isMaybeFacet=FALSE;
-                        P->p1=initialForm->m[jj];       build spolys of initialForm in_v
+                        P->p1=initialForm->m[jj];       //build spolys of initialForm in_v
 
                         for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
                         {
                                 P->p2=initialForm->m[kk];
-                                ksCreateSpoly(P);
-                                if(P->p!=NULL)  spoly non zero=?
-                                {
+                                ksCreateSpoly(P);                                                       
+                                if(P->p!=NULL)  //spoly non zero=?
+                                {       
                                         poly p;//NOTE Don't use pInit here. Evil memleak will follow
                                         poly q;
                                         poly pDel,qDel;
                                         p=pCopy(P->p);
-                                        q=pHead(p);     Monomial q
+                                        q=pHead(p);     //Monomial q
                                         pDelete(&q);
                                         pDel=p; qDel=q;
                                         isMaybeFacet=FALSE;
-                                        TODO: Suffices to check LTs here
+                                        //TODO: Suffices to check LTs here
                                         while(p!=NULL)
-                                        {
+                                        {                                               
                                                 q=pHead(p);
                                                 for(int ll=0;ll<IDELEMS(L);ll++)
                                                 {
                                                         if(pLmEqual(L->m[ll],q) || pDivisibleBy(L->m[ll],q))
-                                                        {
+                                                        {                                                       
                                                                 isMaybeFacet=TRUE;
-                                                                break;for
+                                                                break;//for
                                                         }
                                                 }
@@ -961 +949 @@
                                                 if(isMaybeFacet==TRUE)
                                                 {
-                                                        break;while(p!=NULL)
+                                                        break;//while(p!=NULL)
                                                 }
                                                 p=pNext(p);
-                                        }while
-                                        pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
+                                        }//while
+//                                      pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
                                         if(q!=NULL) pDelete(&q);
                                         pDelete(&pDel);
@@ -971 +959 @@
                                         if(isMaybeFacet==FALSE)
                                         {
-                                                dd_set_si(ddineq->matrix[ii][0],1);
-                                                if(num_alloc==0)
-                                                        num_alloc += 1;
-                                                else
-                                                        num_alloc += 1;
+                                                dd_set_si(ddineq->matrix[ii][0],1);                                             
+//                                              if(num_alloc==0)
+//                                                      num_alloc += 1;
+//                                              else                                            
+//                                                      num_alloc += 1;
                                                 if(num_alloc==num_elts) num_alloc==0 ? num_alloc=1 : num_alloc*=2;
-
+                                                
                                                 void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
                                                 if(!tmp)
@@ -987 +975 @@
                                                 redRowsArray[num_elts]=ii;
                                                 num_elts++;
-                                                break;//for(int kk, since we have found one that is not in L
-                                                goto _start;    mea culpa, mea culpa, mea maxima culpa
+                                                //break;//for(int kk, since we have found one that is not in L  
+                                                goto _start;    //mea culpa, mea culpa, mea maxima culpa
                                         }
-                                }if(P->p!=NULL)
+                                }//if(P->p!=NULL)
                                 pDelete(&(P->p));
-                        }for k
-                }for jj
+                        }//for k
+                }//for jj
                 _start:;
                 idDelete(&L);
                 delete P;
                 idDelete(&initialForm);
-        }for(ii<ddineq-rowsize
-        delete gamma;
-        int offset=0;needed for correction of redRowsArray[ii]
+        }//for(ii<ddineq-rowsize
+//      delete gamma;
+        int offset=0;//needed for correction of redRowsArray[ii]
 #ifdef gfan_DEBUG
         printf("Removed %i of %i in preprocessing step\n",num_elts,ddineq->rowsize);
@@ -1006 +994 @@
         for( int ii=0;ii<num_elts;ii++ )
         {
-                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);cddlib sucks at enumeration
+                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
                 offset++;
-        }
-        free(redRowsArray);NOTE May crash on some machines.
+        }       
+        free(redRowsArray);//NOTE May crash on some machines.
         /*And now for the strictly positive rows
         * Doesn't gain significant speedup
@@ -1022 +1010 @@
                         (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
                 bool isStrictlyPos=FALSE;
-                int posCtr=0;
+                int posCtr=0;           
                 for(int jj=0;jj<this->numVars;jj++)
                 {
@@ -1030 +1018 @@
                         {
                                 if ((*ivPos)[jj]>=0)
-                                {
-                                        posCtr++;
+                                {                               
+                                        posCtr++;                               
                                 }
-                        }
+                        }                       
                         delete ivCanonical;
                 }
@@ -1052 +1040 @@
                         posRowsArray = (int*)tmp;
                         posRowsArray[num_elts]=ii;
-                        num_elts++;
+                        num_elts++;     
                 }
                 delete ivPos;
@@ -1065 +1053 @@
 #endif
 
-        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
-        ddrows = ddineq->rowsize;       Size of the matrix with redundancies removed
+        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);       
+        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
         ddcols = ddineq->colsize;
-
+        
         this->ddFacets = dd_CopyMatrix(ddineq);
-
-        /*Write the normals into class facet*/
-        facet *fAct;    pointer to active facet
+                        
+        /*Write the normals into class facet*/  
+        facet *fAct;    //pointer to active facet       
         int numNonFlip=0;
         for (int kk = 0; kk<ddrows; kk++)
         {
-                int64 ggT=1;NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
-                int64vec *load = new int64vec(this->numVars);int64vec to store a single facet normal that will then be stored via setFacetNormal
+                int64 ggT=1;//NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
+                int64vec *load = new int64vec(this->numVars);//int64vec to store a single facet normal that will then be stored via setFacetNormal
                 for (int jj = 1; jj <ddcols; jj++)
                 {
                         int64 val;
                         val = (int64)mpq_get_d(ddineq->matrix[kk][jj]);
-                        (*load)[jj-1] = val;    store typecasted entry at pos jj-1 of load
+                        (*load)[jj-1] = val;    //store typecasted entry at pos jj-1 of load
                         ggT = int64gcd(ggT,/*(int64&)foo*/val);
-                }for (int jj = 1; jj <ddcols; jj++)
+                }//for (int jj = 1; jj <ddcols; jj++)
                 if(ggT>1)
                 {
                         for(int ll=0;ll<this->numVars;ll++)
-                                (*load)[ll] /= ggT;make primitive vector
+                                (*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
-                */
+                */      
                 if(gcone::hasHomInput==FALSE)
                 {
-                        TODO: No dP needed
+                        //TODO: No dP needed
                         bool isFlip=FALSE;
                         for(int jj = 0; jj<load->length(); jj++)
                         {
-                                int64vec *ivCanonical = new int64vec(load->length());
-                                (*ivCanonical)[jj]=1;
-                                if (dotProduct(*load,*ivCanonical)<0)
-                                {
-                                        isFlip=TRUE;
-                                        break;  //URGHS
-                                }
-                                delete ivCanonical;
+//                              int64vec *ivCanonical = new int64vec(load->length());
+//                              (*ivCanonical)[jj]=1;
+//                              if (dotProduct(*load,*ivCanonical)<0)   
+//                              {
+//                                      isFlip=TRUE;
+//                                      break;  //URGHS
+//                              }
+//                              delete ivCanonical;
                                 if((*load)[jj]<0)
                                 {
                                         isFlip=TRUE;
                                         break;
-                                }
+                                }                               
                         }/*End of check for flippability*/
-                        if(iv64isStrictlyPositive(load))
-                                isFlip=TRUE;
+//                      if(iv64isStrictlyPositive(load))
+//                              isFlip=TRUE;
                         if(isFlip==FALSE)
                         {
@@ -1123 +1111 @@
                                         facet *fRoot = new facet();
                                         this->facetPtr = fRoot;
-                                        fAct = fRoot;
+                                        fAct = fRoot;                                                   
                                 }
                                 else
@@ -1134 +1122 @@
                                 fAct->setUCN(this->getUCN());
 #ifdef gfan_DEBUG
-                                printf("Marking facet (");load->show(1,0);printf(") as non flippable\n");
+                                printf("Marking facet (");load->show(1,0);printf(") as non flippable\n");               
 #endif
                         }
@@ -1153 +1141 @@
                                 fAct->isFlippable=TRUE;
                                 fAct->setFacetNormal(load);
-                                fAct->setUCN(this->getUCN());
+                                fAct->setUCN(this->getUCN());                                   
                         }
-                }hasHomInput==FALSE
-                else    Every facet is flippable
-                {       /*Now load should be full and we can call setFacetNormal*/
+                }//hasHomInput==FALSE
+                else    //Every facet is flippable
+                {       /*Now load should be full and we can call setFacetNormal*/                                      
                         this->numFacets++;
                         if(this->numFacets==1)
@@ -1172 +1160 @@
                         fAct->isFlippable=TRUE;
                         fAct->setFacetNormal(load);
-                        fAct->setUCN(this->getUCN());
-                }if (isFlippable==FALSE)
-                delete load;
-        }for (int kk = 0; kk<ddrows; kk++)
-
-        In cases like I=<x-1,y-1> there are only non-flippable facets...
+                        fAct->setUCN(this->getUCN());                                   
+                }//if (isFlippable==FALSE)
+                delete load;                            
+        }//for (int kk = 0; kk<ddrows; kk++)
+                        
+        //In cases like I=<x-1,y-1> there are only non-flippable facets...
         if(numNonFlip==this->numFacets)
-        {
+        {                                       
                 WerrorS ("Only non-flippable facets. Terminating...\n");
-                exit(-1);//Bit harsh maybe...
-        }
-
+//              exit(-1);//Bit harsh maybe...
+        }
+                        
         /*
         Now we should have a linked list containing the facet normals of those facets that are
@@ -1189 +1177 @@
         -flipable
         Adressing is done via *facetPtr
-        */
+        */                      
         if (compIntPoint==TRUE)
         {
@@ -1198 +1186 @@
                 {
                         dd_set_si(posRestr->matrix[ii][jj],1);
-                        jj++;
+                        jj++;                                                   
                 }
                 dd_MatrixAppendTo(&ddineq,posRestr);
-                interiorPoint(ddineq, *iv);     NOTE ddineq contains non-flippable facets
-                this->setIntPoint(iv);  stores the interior point in gcone::ivIntPt
+                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
+                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
                 delete iv;
                 dd_FreeMatrix(posRestr);
-        }
-        //Clean up but don't delete the return value!
-        dd_FreeMatrix(ddineq);
-        set_free(ddredrows);check
-        set_free(ddlinset);check
-        free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
+        }       
+        //Clean up but don't delete the return value!   
+        //dd_FreeMatrix(ddineq);        
+        set_free(ddredrows);//check
+        set_free(ddlinset);//check
+        //free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
 #ifdef gfanp
         gettimeofday(&end, 0);
@@ -1216 +1204 @@
 #endif
 
-}gcone::getConeNormals(ideal I)
-
+}//gcone::getConeNormals(ideal I)
+                
 /** \brief Compute the (codim-2)-facets of a given cone
  * This method is used during noRevS
@@ -1229 +1217 @@
         gettimeofday(&start, 0);
 #endif
-        this->facetPtr->codim2Ptr = new facet(2);       //instantiate a (codim-2)-facet
+        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
         facet *fAct;
-        fAct = this->facetPtr;
+        fAct = this->facetPtr;          
         facet *codim2Act;
-        codim2Act = this->facetPtr->codim2Ptr;
-        dd_MatrixPtr ddineq;,P,ddakt;
+        //codim2Act = this->facetPtr->codim2Ptr;
+        dd_MatrixPtr ddineq;//,P,ddakt;
         dd_ErrorType err;
-        ddineq = facets2Matrix(gc);     //get a matrix representation of the cone
+        //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
         ddineq = dd_CopyMatrix(gc.ddFacets);
         /*Now set appropriate linearity*/
-        for (int ii=0; ii<this->numFacets; ii++)
-        {
+        for (int ii=0; ii<this->numFacets; ii++)                        
+        {       
                 dd_rowset impl_linset, redset;
                 dd_rowindex newpos;
                 dd_MatrixPtr ddakt;
                 ddakt = dd_CopyMatrix(ddineq);
-                ddakt->representation=dd_Inequality;    //Not using this makes it faster. But why does the quick check below still work?
-                ddakt->representation=dd_Generator;
+//              ddakt->representation=dd_Inequality;    //Not using this makes it faster. But why does the quick check below still work?
+//              ddakt->representation=dd_Generator;
                 set_addelem(ddakt->linset,ii+1);/*Now set appropriate linearity*/
 #ifdef gfanp
                 timeval t_ddMC_start, t_ddMC_end;
                 gettimeofday(&t_ddMC_start,0);
-#endif
-                dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
+#endif                          
+                //dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
                 dd_PolyhedraPtr ddpolyh;
                 ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
-                ddpolyh=dd_DDMatrix2Poly2(ddakt, dd_MaxCutoff, &err);
+//              ddpolyh=dd_DDMatrix2Poly2(ddakt, dd_MaxCutoff, &err);
                 dd_MatrixPtr P;
-                P=dd_CopyGenerators(ddpolyh);
+                P=dd_CopyGenerators(ddpolyh);           
                 dd_FreePolyhedra(ddpolyh);
-                TODO Call for one cone , normalize - check equalities - plus lineality -done
+                //TODO Call for one cone , normalize - check equalities - plus lineality -done
 #ifdef gfanp
                 gettimeofday(&t_ddMC_end,0);
                 t_ddMC += (t_ddMC_end.tv_sec - t_ddMC_start.tv_sec + 1e-6*(t_ddMC_end.tv_usec - t_ddMC_start.tv_usec));
-#endif
+#endif  
                 /* We loop through each row of P normalize it by making all
                 * entries integer ones and add the resulting vector to the
                 * int matrix facet::codim2Facets */
                 for (int jj=1;jj<=/*ddakt*/P->rowsize;jj++)
-                {
+                {                                       
                         fAct->numCodim2Facets++;
-                        if(fAct->numCodim2Facets==1)
-                        {
-                                fAct->codim2Ptr = new facet(2);
+                        if(fAct->numCodim2Facets==1)                                    
+                        {                                               
+                                fAct->codim2Ptr = new facet(2);                                         
                                 codim2Act = fAct->codim2Ptr;
                         }
@@ -1297 +1285 @@
 #endif
                         codim2Act->setFacetNormal(n);
-                        delete n;
-                }
+                        delete n;                                       
+                }               
                 /*We check whether the facet spanned by the codim-2 facets
                 * intersects with the positive orthant. Otherwise we define this
-                * facet to be non-flippable. Works since we set the appropriate
+                * facet to be non-flippable. Works since we set the appropriate 
                 * linearity for ddakt above.
                 */
-                TODO It might be faster to compute jus the implied equations instead of a relative interior point
-                int64vec *iv_intPoint = new int64vec(this->numVars);
-                dd_MatrixPtr shiftMatrix;
-                dd_MatrixPtr intPointMatrix;
-                shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
-                for(int kk=0;kk<this->numVars;kk++)
-                {
-                        dd_set_si(shiftMatrix->matrix[kk][0],1);
-                        dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
-                }
-                intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
- #ifdef gfanp
-                timeval t_iP_start, t_iP_end;
-                gettimeofday(&t_iP_start, 0);
- #endif
-                interiorPoint(intPointMatrix,*iv_intPoint);
- //             dd_rowset impl_linste,lbasis;
- //             dd_LPSolutionPtr lps=NULL;
- //             dd_ErrorType err;
- //             dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
- #ifdef gfanp
-                gettimeofday(&t_iP_end, 0);
-                t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
- #endif
-                for(int ll=0;ll<this->numVars;ll++)
-                {
-                        if( (*iv_intPoint)[ll] < 0 )
-                        {
-                                fAct->isFlippable=FALSE;
-                                break;
-                        }
-                }
+                //TODO It might be faster to compute jus the implied equations instead of a relative interior point
+//              int64vec *iv_intPoint = new int64vec(this->numVars);
+//              dd_MatrixPtr shiftMatrix;
+//              dd_MatrixPtr intPointMatrix;
+//              shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
+//              for(int kk=0;kk<this->numVars;kk++)
+//              {
+//                      dd_set_si(shiftMatrix->matrix[kk][0],1);
+//                      dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
+//              }
+//              intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
+// #ifdef gfanp
+//              timeval t_iP_start, t_iP_end;
+//              gettimeofday(&t_iP_start, 0);
+// #endif               
+//              interiorPoint(intPointMatrix,*iv_intPoint);
+// //           dd_rowset impl_linste,lbasis;
+// //           dd_LPSolutionPtr lps=NULL;
+// //           dd_ErrorType err;
+// //           dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
+// #ifdef gfanp
+//              gettimeofday(&t_iP_end, 0);
+//              t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
+// #endif
+//              for(int ll=0;ll<this->numVars;ll++)
+//              {
+//                      if( (*iv_intPoint)[ll] < 0 )
+//                      {
+//                              fAct->isFlippable=FALSE;
+//                              break;
+//                      }
+//              }
                 /*End of check*/
                 /*This test should be way less time consuming*/
@@ -1358 +1346 @@
                 }
                 if(containsStrictlyPosRay==FALSE)
-                TODO Not sufficient. Intersect with pos orthant for pos int
+                //TODO Not sufficient. Intersect with pos orthant for pos int
                         fAct->isFlippable=FALSE;
 #ifdef gfanp
@@ -1365 +1353 @@
 #endif
                 /**/
-                fAct = fAct->next;
+                fAct = fAct->next;      
                 dd_FreeMatrix(ddakt);
                 dd_FreeMatrix(P);
-        }for
+        }//for  
         dd_FreeMatrix(ddineq);
 #ifdef gfanp
@@ -1375 +1363 @@
 #endif
 }
-
+                
 /** Really extremal rays this time ;)
 * Extremal rays are unique modulo the homogeneity space.
@@ -1383 +1371 @@
 * checking whether the inner product of the ray with the normal is zero.
 * We use ivAdd here which returns a new int64vec. Therefore we need to avoid
-* a memory leak which would be cause by the line
+* a memory leak which would be cause by the line 
 * iv=ivAdd(iv,b)
-* So we keep pointer tmp to iv and delete(tmp), so there should not occur a
+* So we keep pointer tmp to iv and delete(tmp), so there should not occur a 
 * memleak
 * TODO normalization
@@ -1397 +1385 @@
         gettimeofday(&poly_start,0);
 #endif
-        Add lineality space - dd_LinealitySpace
+        //Add lineality space - dd_LinealitySpace
         dd_MatrixPtr ddineq;
-        dd_ErrorType err;
-        if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
-                ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
-        else
-                ddineq = dd_CopyMatrix(gc.ddFacets);
+        dd_ErrorType err;       
+//      if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
+//              ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
+//      else
+//              ddineq = dd_CopyMatrix(gc.ddFacets);
         ddineq = (dd_LinealitySpace->rowsize>0) ? dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace) : dd_CopyMatrix(gc.ddFacets);
         /* In case the input is non-homogeneous we add constrains for the positive orthant.
-        * This is justified by the fact that for non-homog ideals we only consider the
+        * This is justified by the fact that for non-homog ideals we only consider the 
         * restricted fan. This way we can be sure to find strictly positive interior points.
         * This in turn makes life easy when checking for flippability!
@@ -1421 +1409 @@
                 assert(ddineq);
                 dd_FreeMatrix(ddPosRestr);
-        }
+        }       
         dd_PolyhedraPtr ddPolyh;
         ddPolyh = dd_DDMatrix2Poly(ddineq, &err);
@@ -1434 +1422 @@
         /* Compute interior point on the fly*/
         int64vec *ivIntPointOfCone = new int64vec(this->numVars);
-        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++)
-        {
-                mpq_init(colSum[jj]);
-                for(int ii=0;ii<P->rowsize;ii++)
-                {
-                        mpq_t tmp; mpq_init(tmp);
-                        mpq_t sum; mpq_init(sum);
-                        mpq_set(sum,colSum[jj]);
-                        mpq_add(tmp,sum,P->matrix[ii][jj+1]);
-                        mpq_set(colSum[jj],tmp);
-                        mpq_clear(tmp);
-                        mpq_clear(sum);
-                }
-                mpz_t den; mpz_init(den);
-                mpq_get_den(den,colSum[jj]);
-                denom[jj]=(int)mpz_get_si(den);
-                mpz_clear(den);
-        }
-        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);
-        mpz_t den; mpz_init(den);
-        mpz_t tmp; mpz_init(tmp);
-        mpq_get_den(tmp,colSum[0]);
-        for (int ii=0;ii<(this->numVars)-1;ii++)
-        {
-                mpq_get_den(den,colSum[ii+1]);
-                mpz_lcm(kgV,tmp,den);
-                mpz_set(tmp, kgV);
-        }
-        mpq_t qkgV;
-        mpq_init(qkgV);
-        mpq_set_z(qkgV,kgV);
-        mpz_clear(kgV);
-        mpz_clear(den);
-        mpz_clear(tmp);*/
         int64vec *foo = new int64vec(this->numVars);
         for(int ii=0;ii<P->rowsize;ii++)
         {
-                int64vec *foo = new int64vec(this->numVars);
+//              int64vec *foo = new int64vec(this->numVars);
                 int64vec *tmp = ivIntPointOfCone;
                 makeInt(P,ii+1,*foo);
                 ivIntPointOfCone = iv64Add(ivIntPointOfCone,foo);
                 delete tmp;
-                delete foo;
+//              delete foo;
         }
         delete foo;
@@ -1488 +1436 @@
         for (int ii=0;ii<(this->numVars);ii++)
         {
-                mpq_t product;
-                mpq_init(product);
-                mpq_mul(product,qkgV,colSum[ii]);
-                (*ivIntPointOfCone)[ii]=(int64)mpz_get_d(mpq_numref(product));
-                if( (*ivIntPointOfCone)[ii]>INT_MAX )
+//              mpq_t product;
+//              mpq_init(product);
+//              mpq_mul(product,qkgV,colSum[ii]);
+//              (*ivIntPointOfCone)[ii]=(int64)mpz_get_d(mpq_numref(product));
+                if( (*ivIntPointOfCone)[ii]>INT_MAX ) 
                         WarnS("Interior point exceeds INT_MAX!\n");
-                mpq_clear(product);
-                Compute intgcd
+//              mpq_clear(product);
+                //Compute intgcd
                 ggT=int64gcd(ggT,(*ivIntPointOfCone)[ii]);
         }
-
-        Divide out a common gcd > 1
+        
+        //Divide out a common gcd > 1
         if(ggT>1)
         {
@@ -1508 +1456 @@
                 }
         }
-        mpq_clear(qkgV);
-        delete [] colSum;
+//      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)
@@ -1517 +1465 @@
                 for(int ii=0;ii<this->numVars;ii++)
                 {
-                        (*ivOne)[ii]=1;
+//                      (*ivOne)[ii]=1;
                         if ((*ivIntPointOfCone)[ii]<maxNegEntry) maxNegEntry=(*ivIntPointOfCone)[ii];
                 }
                 maxNegEntry *= -1;
-                maxNegEntry++;To be on the safe side
+                maxNegEntry++;//To be on the safe side
                 for(int ii=0;ii<this->numVars;ii++)
                         (*ivOne)[ii]=maxNegEntry;
@@ -1527 +1475 @@
                 ivIntPointOfCone=iv64Add(ivIntPointOfCone,ivOne);
                 delete(tmp);
-                while( !iv64isStrictlyPositive(ivIntPointOfCone) )
-                {
-                        int64vec *tmp = ivIntPointOfCone;
-                        for(int jj=0;jj<this->numVars;jj++)
-                                (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
-                        ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
-                        delete tmp;
-                }
+//              while( !iv64isStrictlyPositive(ivIntPointOfCone) )
+//              {
+//                      int64vec *tmp = ivIntPointOfCone;
+//                      for(int jj=0;jj<this->numVars;jj++)
+//                              (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
+//                      ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
+//                      delete tmp;                             
+//              }
                 delete ivOne;
                 int64 ggT=(*ivIntPointOfCone)[0];
@@ -1545 +1493 @@
                 }
         }
-        assert(iv64isStrictlyPositive(ivIntPointOfCone));
-
+//      assert(iv64isStrictlyPositive(ivIntPointOfCone));
+        
         this->setIntPoint(ivIntPointOfCone);
         delete(ivIntPointOfCone);
         /* end of interior point computation*/
-
-        Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
+        
+        //Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
         int rows=P->rowsize;
         facet *fAct=gc.facetPtr;
-        Construct an array to hold the extremal rays of the cone
-        this->gcRays = (int64vec**)omAlloc0(sizeof(int64vec*)*P->rowsize);
+        //Construct an array to hold the extremal rays of the cone
+        this->gcRays = (int64vec**)omAlloc0(sizeof(int64vec*)*P->rowsize);      
         for(int ii=0;ii<P->rowsize;ii++)
         {
                 int64vec *rowvec = new int64vec(this->numVars);
-                makeInt(P,ii+1,*rowvec);get an integer entry instead of rational, rowvec is primitve
+                makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve
                 this->gcRays[ii] = iv64Copy(rowvec);
                 delete rowvec;
-        }
+        }       
         this->numRays=P->rowsize;
-        Check which rays belong to which facet
+        //Check which rays belong to which facet
         while(fAct!=NULL)
         {
-                const int64vec *fNormal; = new int64vec(this->numVars);
-                fNormal = fAct->getRef2FacetNormal();->getFacetNormal();
+                const int64vec *fNormal;// = new int64vec(this->numVars);
+                fNormal = fAct->getRef2FacetNormal();//->getFacetNormal();
                 int64vec *ivIntPointOfFacet = new int64vec(this->numVars);
                 for(int ii=0;ii<rows;ii++)
-                {
+                {                       
                         if(dotProduct(*fNormal,this->gcRays[ii])==0)
                         {
-                                int64vec *tmp = ivIntPointOfFacet;Prevent memleak
+                                int64vec *tmp = ivIntPointOfFacet;//Prevent memleak
                                 fAct->numCodim2Facets++;
                                 facet *codim2Act;
-                                if(fAct->numCodim2Facets==1)
-                                {
-                                        fAct->codim2Ptr = new facet(2);
+                                if(fAct->numCodim2Facets==1)                                    
+                                {                                               
+                                        fAct->codim2Ptr = new facet(2);                                         
                                         codim2Act = fAct->codim2Ptr;
                                 }
@@ -1587 +1535 @@
                                         codim2Act = codim2Act->next;
                                 }
-                                codim2Act->setFacetNormal(rowvec);
-                                Rather just let codim2Act point to the corresponding int64vec of gcRays
+                                //codim2Act->setFacetNormal(rowvec);
+                                //Rather just let codim2Act point to the corresponding int64vec of gcRays
                                 codim2Act->fNormal=this->gcRays[ii];
-                                fAct->numRays++;
-                                Memleak avoided via tmp
+                                fAct->numRays++;                                 
+                                //Memleak avoided via tmp
                                 ivIntPointOfFacet=iv64Add(ivIntPointOfFacet,this->gcRays[ii]);
-                                Now tmp still points to the OLD address of ivIntPointOfFacet
+                                //Now tmp still points to the OLD address of ivIntPointOfFacet
                                 delete(tmp);
-
+                                        
                         }
-                }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)
+                }//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)
                 if( !iv64isStrictlyPositive(ivIntPointOfFacet) && hasHomInput==TRUE)
                 {
@@ -1611 +1559 @@
                                 for(int jj=0;jj<this->numVars;jj++)
                                 {
-                                        (*ivOne)[jj] = (*ivOne)[jj] << 1; times 2
+                                        (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
                                 }
                                 ivIntPointOfFacet = iv64Add(ivIntPointOfFacet/*diff*/,ivOne);
-                                delete tmp;
+                                delete tmp;                             
                         }
                         delete ivOne;
@@ -1625 +1573 @@
                         for(int ii=0;ii<this->numVars;ii++)
                                 (*ivIntPointOfFacet)[ii] /= ggT;
-                }
+                }                       
                 fAct->setInteriorPoint(ivIntPointOfFacet);
-
+                
                 delete(ivIntPointOfFacet);
-                Now (if we have at least 3 variables) do a bubblesort on the rays
+                //Now (if we have at least 3 variables) do a bubblesort on the rays
                 /*if(this->numVars>2)
                 {
@@ -1643 +1591 @@
                         do
                         {
-                                exchanged=FALSE;n=fAct->numRays-1;
+                                exchanged=FALSE;//n=fAct->numRays-1;                            
                                 for(unsigned ii=0;ii<=n-1;ii++)
-                                {
+                                {                                       
                                         if((A[ii]->fNormal)->compare((A[ii+1]->fNormal))==1)
                                         {
-                                                Swap rays
+                                                //Swap rays
                                                 cout << "Swapping ";
                                                 A[ii]->fNormal->show(1,0); cout << " with "; A[ii+1]->fNormal->show(1,0); cout << endl;
@@ -1657 +1605 @@
                                                 if(ii==0)
                                                         fAct->codim2Ptr=A[ii+1];
-                                                end swap
-                                                facet *tmp=A[ii];swap in list
+                                                //end swap
+                                                facet *tmp=A[ii];//swap in list
                                                 A[ii+1]=A[ii];
                                                 A[ii]=tmp;
-                                                tmp=NULL;
-                                        }
+//                                              tmp=NULL;                       
+                                        }                                       
                                 }
-                                n--;
+                                n--;                    
                         }while(exchanged==TRUE && n>=0);
-                }*/if pVariables>2
-                delete fNormal;
+                }*///if pVariables>2
+//              delete fNormal;         
                 fAct = fAct->next;
-        }end of facet checking
+        }//end of facet checking
         dd_FreeMatrix(P);
-        Now all extremal rays should be set w.r.t their respective fNormal
-        TODO Not sufficient -> vol2 II/125&127
-        NOTE Sufficient according to cddlibs doc. These ARE rays
+        //Now all extremal rays should be set w.r.t their respective fNormal
+        //TODO Not sufficient -> vol2 II/125&127
+        //NOTE Sufficient according to cddlibs doc. These ARE rays
         //What the hell... let's just take interior points
         if(gcone::hasHomInput==FALSE)
@@ -1680 +1628 @@
                 while(fAct!=NULL)
                 {
-                        bool containsStrictlyPosRay=FALSE;
-                        facet *codim2Act;
-                        codim2Act = fAct->codim2Ptr;
-                        while(codim2Act!=NULL)
-                        {
-                                int64vec *rayvec;
-                                rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
-                                //int negCtr=0;
-                                if(iv64isStrictlyPositive(rayvec))
-                                {
-                                        containsStrictlyPosRay=TRUE;
-                                        delete(rayvec);
-                                        break;
-                                }
-                                delete(rayvec);
-                                codim2Act = codim2Act->next;
-                        }
-                        if(containsStrictlyPosRay==FALSE)
-                                fAct->isFlippable=FALSE;
+//                      bool containsStrictlyPosRay=FALSE;
+//                      facet *codim2Act;
+//                      codim2Act = fAct->codim2Ptr;
+//                      while(codim2Act!=NULL)
+//                      {
+//                              int64vec *rayvec;
+//                              rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
+//                              //int negCtr=0;
+//                              if(iv64isStrictlyPositive(rayvec))
+//                              {
+//                                      containsStrictlyPosRay=TRUE;
+//                                      delete(rayvec);
+//                                      break;
+//                              }                               
+//                              delete(rayvec);
+//                              codim2Act = codim2Act->next;
+//                      }
+//                      if(containsStrictlyPosRay==FALSE)
+//                              fAct->isFlippable=FALSE;
                         if(!iv64isStrictlyPositive(fAct->interiorPoint))
                                 fAct->isFlippable=FALSE;
                         fAct = fAct->next;
                 }
-        }hasHomInput?
+        }//hasHomInput? 
 #ifdef gfanp
         gettimeofday(&end, 0);
         t_getExtremalRays += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
-#endif
+#endif  
 }
 
@@ -1713 +1661 @@
 void gcone::orderRays()
 {
-   qsort(gcRays,sizeof(int64vec),int64vec::compare);
-}
-
+//   qsort(gcRays,sizeof(int64vec),int64vec::compare);
+}
+        
 inline bool gcone::iv64isStrictlyPositive(const int64vec * iv64)
 {
@@ -1725 +1673 @@
                         res=FALSE;
                         break;
-                }
+                }               
         }
         return res;
 }
 
-/** \brief Compute the Groebner Basis on the other side of a shared facet
+/** \brief Compute the Groebner Basis on the other side of a shared facet 
  *
  * Implements algorithm 4.3.2 from Anders' thesis.
@@ -1738 +1686 @@
  * Parallelity is checked using basic linear algebra. See gcone::isParallel.
  * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
- * computing an interior point of the facet and taking all terms having the same weight with respect
+ * computing an interior point of the facet and taking all terms having the same weight with respect 
  * to this interior point.
  *\param ideal, facet
  * Input: a marked,reduced Groebner basis and a facet
  */
-inline void gcone::flip(ideal gb, facet *f)             Compute "the other side"
-{
+inline void gcone::flip(ideal gb, facet *f)             //Compute "the other side"
+{       
 #ifdef gfanp
         timeval start, end;
         gettimeofday(&start, 0);
-#endif
-        int64vec *fNormal; = new int64vec(this->numVars);       //facet normal, check for parallelity
-        fNormal = f->getFacetNormal();  read this->fNormal;
+#endif          
+        int64vec *fNormal;// = new int64vec(this->numVars);     //facet normal, check for parallelity                   
+        fNormal = f->getFacetNormal();  //read this->fNormal;
 #ifdef gfan_DEBUG
-        std::cout << "running gcone::flip" << std::endl;
+//      std::cout << "running gcone::flip" << std::endl;
         printf("flipping UCN %i over facet",this->getUCN());
         fNormal->show(1,0);
-        printf(") with UCN %i\n",f->getUCN() );
+        printf(") with UCN %i\n",f->getUCN() ); 
 #endif
         if(this->getUCN() != f->getUCN())
@@ -1763 +1711 @@
         }
         /*1st step: Compute the initial ideal*/
-        /*poly initialFormElement[IDELEMS(gb)];*/       array of #polys in GB to store initial form
+        /*poly initialFormElement[IDELEMS(gb)];*/       //array of #polys in GB to store initial form
         ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
-
+        
         computeInv(gb,initialForm,*fNormal);
 
@@ -1771 +1719 @@
 /*      cout << "Initial ideal is: " << endl;
         idShow(initialForm);
-        f->printFlipGB();*/
-        cout << "===" << endl;
-#endif
+        //f->printFlipGB();*/
+//      cout << "===" << endl;
+#endif                  
         /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
         /*Substep 2.1
-        compute $G_{-\alpha}(in_v(I))
+        compute $G_{-\alpha}(in_v(I)) 
         see journal p. 66
         NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
@@ -1783 +1731 @@
         ring srcRing=currRing;
         ring tmpRing;
-
+                        
         if( (srcRing->order[0]!=ringorder_a))
         {
-                int64vec *iv; = new int64vec(this->numVars);
-                iv = ivNeg(fNormal);ivNeg uses iv64Copy -> new
-                tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
+                int64vec *iv;// = new int64vec(this->numVars);
+                iv = ivNeg(fNormal);//ivNeg uses iv64Copy -> new
+//              tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
                 tmpRing=rCopyAndAddWeight(srcRing,iv);
                 delete iv;
@@ -1805 +1753 @@
                 tmpRing->block1[0]=length;
                 rComplete(tmpRing);
-                omFree(A);
-        }
-        delete fNormal;
-        rChangeCurrRing(tmpRing);
-
-        ideal ina;
+                //omFree(A);
+        }
+        delete fNormal; 
+        rChangeCurrRing(tmpRing);       
+                        
+        ideal ina;                      
         ina=idrCopyR(initialForm,srcRing);
         idDelete(&initialForm);
         ideal H;
-        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
+//      H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
 #ifdef gfanp
         timeval t_kStd_start, t_kStd_end;
@@ -1822 +1770 @@
                 H=kStd(ina,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
         else
-                H=kStd(ina,NULL,isNotHomog,NULL);       This is \mathcal(G)_{>_-\alpha}(in_v(I))
+                H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
 #ifdef gfanp
         gettimeofday(&t_kStd_end, 0);
@@ -1835 +1783 @@
         rChangeCurrRing(srcRing);
         ideal srcRing_H;
-        ideal srcRing_HH;
+        ideal srcRing_HH;                       
         srcRing_H=idrCopyR(H,tmpRing);
         //H is needed further below, so don't idDelete here
         srcRing_HH=ffG(srcRing_H,this->gcBasis);
         idDelete(&srcRing_H);
-
+                
         /*Substep 2.2.1
          * Mark according to G_-\alpha
@@ -1846 +1794 @@
          * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
          * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
-         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
+         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we 
          * compute the difference accordingly
         */
@@ -1852 +1800 @@
         timeval t_markings_start, t_markings_end;
         gettimeofday(&t_markings_start, 0);
-#endif
+#endif          
         bool markingsAreCorrect=FALSE;
         dd_MatrixPtr intPointMatrix;
         int iPMatrixRows=0;
-        dd_rowrange aktrow=0;
+        dd_rowrange aktrow=0;                   
         for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
         {
                 poly aktpoly=(poly)srcRing_HH->m[ii];//This is a pointer, so don't pDelete
-                iPMatrixRows = iPMatrixRows+pLength(aktpoly);
+                iPMatrixRows = iPMatrixRows+pLength(aktpoly);           
         }
         /* additionally one row for the standard-simplex and another for a row that becomes 0 during
          * construction of the differences
          */
-        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
-        intPointMatrix->numbtype=dd_Integer;    NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
-
+        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); 
+        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
+        
         for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
         {
-                markingsAreCorrect=FALSE;       crucial to initialise here
+                markingsAreCorrect=FALSE;       //crucial to initialise here
                 poly aktpoly=srcRing_HH->m[ii]; //Only a pointer, so don't pDelete
                 /*Comparison of leading monomials is done via exponent vectors*/
@@ -1878 +1826 @@
                         int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
                         pGetExpV(aktpoly,src_ExpV);
-                        rChangeCurrRing(tmpRing);       this ring change is crucial!
+                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
                         poly p=pCopy(H->m[ii]);
                         pGetExpV(p/*pCopy(H->m[ii])*/,dst_ExpV);
@@ -1887 +1835 @@
                         {
 #ifdef gfan_DEBUG
-                                cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
+//                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
 #endif
                                 if (src_ExpV[kk]!=dst_ExpV[kk])
@@ -1893 +1841 @@
                                         expVAreEqual=FALSE;
                                 }
-                        }
+                        }                                                               
                         if (expVAreEqual==TRUE)
                         {
-                                markingsAreCorrect=TRUE; everything is fine
+                                markingsAreCorrect=TRUE; //everything is fine
 #ifdef gfan_DEBUG
-                                cout << "correct markings" << endl;
-#endif
-                        }if (pHead(aktpoly)==pHead(H->m[jj])
+//                              cout << "correct markings" << endl;
+#endif
+                        }//if (pHead(aktpoly)==pHead(H->m[jj])
                         omFree(src_ExpV);
                         omFree(dst_ExpV);
-                }for (int jj=0;jj<IDELEMS(H);jj++)
-
+                }//for (int jj=0;jj<IDELEMS(H);jj++)
+                
                 int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
                 if (markingsAreCorrect==TRUE)
@@ -1913 +1861 @@
                 {
                         rChangeCurrRing(tmpRing);
-                        pGetExpV(pHead(H->m[ii]),leadExpV); We use H->m[ii] as leading monomial
+                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
                         rChangeCurrRing(srcRing);
                 }
-                /*compute differences of the expvects*/
+                /*compute differences of the expvects*/                         
                 while (pNext(aktpoly)!=NULL)
                 {
                         int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
-                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
+                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term) 
                         is not omitted when computing the differences*/
                         if(markingsAreCorrect==TRUE)
@@ -1934 +1882 @@
                         int ctr=0;
                         for (int jj=0;jj<this->numVars;jj++)
-                        {
-                                /*Store into ddMatrix*/
+                        {                               
+                                /*Store into ddMatrix*/                         
                                 if(leadExpV[jj+1]-v[jj+1])
                                         ctr++;
@@ -1941 +1889 @@
                         }
                         /*It ought to be more sensible to avoid 0-rows in the first place*/
-                        if(ctr==this->numVars)//We have a 0-row
-                                dd_MatrixRowRemove(&intPointMatrix,aktrow);
-                        else
+//                      if(ctr==this->numVars)//We have a 0-row
+//                              dd_MatrixRowRemove(&intPointMatrix,aktrow);
+//                      else
                                 aktrow +=1;
                         omFree(v);
                 }
                 omFree(leadExpV);
-        }for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
+        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
 #ifdef gfanp
         gettimeofday(&t_markings_end, 0);
@@ -1960 +1908 @@
         preprocessInequalities(intPointMatrix);
         /*Now we add the constraint for the standard simplex*/
-        dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
-        for (int jj=1;jj<=this->numVars;jj++)
-        {
-                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
-        }
+//      dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
+//      for (int jj=1;jj<=this->numVars;jj++)
+//      {
+//              dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
+//      }
         //Let's make sure we compute interior points from the positive orthant
-        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
-
-        int jj=1;
-        for (int ii=0;ii<this->numVars;ii++)
-        {
-                dd_set_si(posRestr->matrix[ii][jj],1);
-                jj++;
-        }
+//      dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
+//      
+//      int jj=1;
+//      for (int ii=0;ii<this->numVars;ii++)
+//      {
+//              dd_set_si(posRestr->matrix[ii][jj],1);
+//              jj++;                                                   
+//      }
         /*We create a matrix containing the standard simplex
         * and constraints to assure a strictly positive point
@@ -1983 +1931 @@
                 {
                         dd_set_si(posRestr->matrix[ii][0],-1);
-                        for(int jj=1;jj<=this->numVars;jj++)
-                                dd_set_si(posRestr->matrix[ii][jj],1);
+                        for(int jj=1;jj<=this->numVars;jj++)                    
+                                dd_set_si(posRestr->matrix[ii][jj],1);                  
                 }
                 else
@@ -2005 +1953 @@
         dd_rowindex newpos;
 
-        NOTE Here we should remove interiorPoint and instead
-         create and ordering like (a(omega),a(fNormal),dp)
-        if(this->ivIntPt==NULL)
-                interiorPoint(intPointMatrix, *iv_weight);      iv_weight now contains the interior point
-        else
-                iv_weight=this->getIntPoint();
+        //NOTE Here we should remove interiorPoint and instead
+        // create and ordering like (a(omega),a(fNormal),dp)
+//      if(this->ivIntPt==NULL)
+                interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
+//      else
+//              iv_weight=this->getIntPoint();
         dd_FreeMatrix(intPointMatrix);
         /*Crude attempt for interior point */
@@ -2034 +1982 @@
         gettimeofday(&t_dd_end, 0);
         t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec));
-#endif
-
+#endif  
+        
         /*Step 3
-        * turn the minimal basis into a reduced one */
-         NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
-         Thus:
-        ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
+        * turn the minimal basis into a reduced one */                  
+        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
+        // Thus: 
+        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight); 
         ring dstRing=rCopy0(tmpRing);
         int length=iv_weight->length();
@@ -2054 +2002 @@
         delete iv_weight;
 
-        ideal dstRing_I;
+        ideal dstRing_I;                        
         dstRing_I=idrCopyR(srcRing_HH,srcRing);
-        idDelete(&srcRing_HH); Hmm.... causes trouble - no more
-        dstRing_I=idrCopyR(inputIdeal,srcRing);
+        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more
+        //dstRing_I=idrCopyR(inputIdeal,srcRing);
         BITSET save=test;
         test|=Sy_bit(OPT_REDSB);
         test|=Sy_bit(OPT_REDTAIL);
 #ifdef gfan_DEBUG
-        test|=Sy_bit(6);        //OPT_DEBUG
+//      test|=Sy_bit(6);        //OPT_DEBUG
 #endif
         ideal tmpI;
-        NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
-        tmpI = idrCopyR(this->inputIdeal,this->baseRing);
+        //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
+        //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
         tmpI = dstRing_I;
 #ifdef gfanp
@@ -2080 +2028 @@
 #endif
         idDelete(&tmpI);
-        idNorm(dstRing_I);
-        kInterRed(dstRing_I);
+        idNorm(dstRing_I);                      
+//      kInterRed(dstRing_I);
         idSkipZeroes(dstRing_I);
         test=save;
         /*End of step 3 - reduction*/
-
-        f->setFlipGB(dstRing_I);store the flipped GB
-        idDelete(&dstRing_I);
-        f->flipRing=rCopy(dstRing);     store the ring on the other side
+                        
+        f->setFlipGB(dstRing_I);//store the flipped GB
+//      idDelete(&dstRing_I);
+        f->flipRing=rCopy(dstRing);     //store the ring on the other side
 #ifdef gfan_DEBUG
         printf("Flipped GB is UCN %i:\n",counter+1);
         idDebugPrint(dstRing_I);
         printf("\n");
-#endif
-        idDelete(&dstRing_I);
-        rChangeCurrRing(srcRing);       return to the ring we started the computation of flipGB in
+#endif  
+        idDelete(&dstRing_I);   
+        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
         rDelete(dstRing);
 #ifdef gfanp
@@ -2101 +2049 @@
         time_flip += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
 #endif
-}void flip(ideal gb, facet *f)
+}//void flip(ideal gb, facet *f)
 
 /** \brief A slightly different approach to flipping
@@ -2111 +2059 @@
 * will be from a ring with (a(),dp,C) and our resulting cone from (a(),a(),dp,C). Hence a way
 * must be found to circumvent the sequence of a()'s growing to a ridiculous size.
-* Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we
-* do have an interior point of the cone by adding the extremal rays. So we replace
-* the latter cone by a cone with (a(sum_of_rays),dp,C).
+* Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we 
+* do have an interior point of the cone by adding the extremal rays. So we replace 
+* the latter cone by a cone with (a(sum_of_rays),dp,C). 
 * Con: It's incredibly ugly
 * Pro: No messing around with readConeFromFile()
@@ -2125 +2073 @@
 #endif
         const int64vec *fNormal;
-        fNormal = f->getRef2FacetNormal();/*->getFacetNormal();*/       read this->fNormal;
+        fNormal = f->getRef2FacetNormal();/*->getFacetNormal();*/       //read this->fNormal;
 #ifdef gfan_DEBUG
         printf("flipping UCN %i over facet(",this->getUCN());
         fNormal->show(1,0);
-        printf(") with UCN %i\n",f->getUCN());
+        printf(") with UCN %i\n",f->getUCN());  
 #endif
         if(this->getUCN() != f->getUCN())
@@ -2137 +2085 @@
         }
         /*1st step: Compute the initial ideal*/
-        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
+        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);      
         computeInv( gb, initialForm, *fNormal );
         ring srcRing=currRing;
         ring tmpRing;
-
+        
         const int64vec *intPointOfFacet;
         intPointOfFacet=f->getInteriorPoint();
-        Now we need two blocks of ringorder_a!
-        May assume the same situation as in flip() here
+        //Now we need two blocks of ringorder_a!
+        //May assume the same situation as in flip() here                       
         if( (srcRing->order[0]!=ringorder_a/*64*/) && (srcRing->order[1]!=ringorder_a/*64*/) )
         {
-                int64vec *iv = new int64vec(this->numVars);init with 1s, since we do not need a 2nd block here but later
-                int64vec *iv_foo = new int64vec(this->numVars,1);//placeholder
-                int64vec *ivw = ivNeg(const_cast<int64vec*>(fNormal));
+                int64vec *iv = new int64vec(this->numVars);//init with 1s, since we do not need a 2nd block here but later
+//              int64vec *iv_foo = new int64vec(this->numVars,1);//placeholder
+                int64vec *ivw = ivNeg(const_cast<int64vec*>(fNormal));          
                 tmpRing=rCopyAndAddWeight2(srcRing,ivw/*intPointOfFacet*/,iv);
                 delete iv;delete ivw;
-                delete iv_foo;
-        }
-        else
+//              delete iv_foo;
+        }
+        else 
         {
                 int64vec *iv=new int64vec(this->numVars);
@@ -2168 +2116 @@
                 {
                         A1[jj] = -(*fNormal)[jj];
-                        A2[jj] = 1;-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
+                        A2[jj] = 1;//-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
                 }
                 omFree(tmpRing->wvhdl[0]);
                 if(tmpRing->wvhdl[1]!=NULL)
                         omFree(tmpRing->wvhdl[1]);
-                tmpRing->wvhdl[0]=(int*)A1;
+                tmpRing->wvhdl[0]=(int*)A1;             
                 tmpRing->block1[0]=length;
-                tmpRing->wvhdl[1]=(int*)A2;
+                tmpRing->wvhdl[1]=(int*)A2;             
                 tmpRing->block1[1]=length;
                 rComplete(tmpRing);*/
         }
-        delete fNormal; //NOTE Do not delete when using getRef2FacetNormal();
-        rChangeCurrRing(tmpRing);
-        Now currRing should have (a(),a(),dp,C)
-        ideal ina;
+//      delete fNormal; //NOTE Do not delete when using getRef2FacetNormal();
+        rChangeCurrRing(tmpRing);       
+        //Now currRing should have (a(),a(),dp,C)               
+        ideal ina;                      
         ina=idrCopyR(initialForm,srcRing);
         idDelete(&initialForm);
@@ -2193 +2141 @@
         test|=Sy_bit(OPT_REDSB);
         test|=Sy_bit(OPT_REDTAIL);
-        if(gcone::hasHomInput==TRUE)
+//      if(gcone::hasHomInput==TRUE)
                 H=kStd(ina,NULL,testHomog/*isHomog*/,NULL/*,gcone::hilbertFunction*/);
-        else
-                H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
+//      else
+//              H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
         test=save;
 #ifdef gfanp
@@ -2204 +2152 @@
         idSkipZeroes(H);
         idDelete(&ina);
-
+        
         rChangeCurrRing(srcRing);
         ideal srcRing_H;
-        ideal srcRing_HH;
+        ideal srcRing_HH;                       
         srcRing_H=idrCopyR(H,tmpRing);
         //H is needed further below, so don't idDelete here
         srcRing_HH=ffG(srcRing_H,this->gcBasis);
         idDelete(&srcRing_H);
-        Now BBA(srcRing_HH) with (a(),a(),dp)
+        //Now BBA(srcRing_HH) with (a(),a(),dp)
         /* Evil modification of currRing */
         ring dstRing=rCopy0(tmpRing);
@@ -2222 +2170 @@
         {
                 A1[jj] = (*intPointOfFacet)[jj];
-                A2[jj] = -(*ivw)[jj];TODO Or minus (*ivw)[ii] ??? NOTE minus
+                A2[jj] = -(*ivw)[jj];//TODO Or minus (*ivw)[ii] ??? NOTE minus
         }
         omFree(dstRing->wvhdl[0]);
         if(dstRing->wvhdl[1]!=NULL)
                 omFree(dstRing->wvhdl[1]);
-        dstRing->wvhdl[0]=(int*)A1;
+        dstRing->wvhdl[0]=(int*)A1;             
         dstRing->block1[0]=length;
-        dstRing->wvhdl[1]=(int*)A2;
+        dstRing->wvhdl[1]=(int*)A2;             
         dstRing->block1[1]=length;
         rComplete(dstRing);
         rChangeCurrRing(dstRing);
-        ideal dstRing_I;
+        ideal dstRing_I;                        
         dstRing_I=idrCopyR(srcRing_HH,srcRing);
-        idDelete(&srcRing_HH); Hmm.... causes trouble - no more
+        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more       
         save=test;
         test|=Sy_bit(OPT_REDSB);
@@ -2242 +2190 @@
         tmpI = dstRing_I;
 #ifdef gfanp
-        timeval t_kStd_start, t_kStd_end;
+//      timeval t_kStd_start, t_kStd_end;
         gettimeofday(&t_kStd_start,0);
 #endif
-        if(gcone::hasHomInput==TRUE)
-                dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
-        else
+//      if(gcone::hasHomInput==TRUE)
+//              dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
+//      else
                 dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
 #ifdef gfanp
@@ -2254 +2202 @@
 #endif
         idDelete(&tmpI);
-        idNorm(dstRing_I);
+        idNorm(dstRing_I);                      
         idSkipZeroes(dstRing_I);
         test=save;
         /*End of step 3 - reduction*/
-
+        
         f->setFlipGB(dstRing_I);
         f->flipRing=rCopy(dstRing);
         rDelete(tmpRing);
         rDelete(dstRing);
-        Now we should have dstRing with (a(),a(),dp,C)
-        This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
-        of cones in noRevS
+        //Now we should have dstRing with (a(),a(),dp,C)
+        //This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
+        //of cones in noRevS
         rChangeCurrRing(srcRing);
 #ifdef gfanp
@@ -2271 +2219 @@
         time_flip2 += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
 #endif
-}flip2
+}//flip2
 
 /** \brief Compute initial ideal
@@ -2289 +2237 @@
                 poly aktpoly = (poly)gb->m[ii];//Ptr, so don't pDelete(aktpoly)
                 int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
-                pGetExpV(aktpoly,leadExpV);     find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
+                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
                 initialFormElement=pHead(aktpoly);
-                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
+//              int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
                 while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
                 {
                         int64vec *check = new int64vec(this->numVars);
-                        aktpoly=pNext(aktpoly); next term
+                        aktpoly=pNext(aktpoly); //next term
                         int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
-                        pGetExpV(aktpoly,v);
+                        pGetExpV(aktpoly,v);            
                         /* Convert (int)v into (int64vec)check */
-                        bool notPar=FALSE;
+//                      bool notPar=FALSE;
                         for (int jj=0;jj<this->numVars;jj++)
                         {
                                 (*check)[jj]=v[jj+1]-leadExpV[jj+1];
-                                register int64 foo=(fNormal)[jj];
-                                if( ( (*check)[jj]  == /*fNormal[jj]*/foo )
-                                 || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) )
-                                {
-                                        notPar=TRUE;
-                                        break;
-                                }
+//                              register int64 foo=(fNormal)[jj];
+//                              if( ( (*check)[jj]  == /*fNormal[jj]*/foo ) 
+//                               || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) ) 
+//                              {
+//                                      notPar=TRUE;
+//                                      break;
+//                              }
                         }
                         omFree(v);
-                        if (isParallel(*check,fNormal))Found a parallel vector. Add it
-                        if(notPar==FALSE)
-                        {
-                                initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));pAdd = p_Add_q destroys args
+                        if (isParallel(*check,fNormal))//Found a parallel vector. Add it
+//                      if(notPar==FALSE)
+                        {
+                                initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));//pAdd = p_Add_q destroys args
                         }
                         delete check;
-                }while
-                omFree(v);
+                }//while
+//              omFree(v);
 #ifdef gfan_DEBUG
-                cout << "Initial Form=";
-                pWrite(initialFormElement[ii]);
-                cout << "---" << endl;
+//              cout << "Initial Form=";                                
+//              pWrite(initialFormElement[ii]);
+//              cout << "---" << endl;
 #endif
                 /*Now initialFormElement must be added to (ideal)initialForm */
@@ -2329 +2277 @@
                 pDelete(&initialFormElement);
                 omFree(leadExpV);
-        }for
+        }//for
 #ifdef gfanp
         gettimeofday(&end, 0);
@@ -2343 +2291 @@
  * compute the factors \f$ a_i \f$
  */
-NOTE: Should be replaced by kNF or kNF2
-NOTE: Done
-NOTE: removed with r12286
-
+//NOTE: Should be replaced by kNF or kNF2
+//NOTE: Done
+//NOTE: removed with r12286
+                
 /** \brief Compute \f$ f-f^{\mathcal{G}} \f$
 */
-NOTE: use kNF or kNF2 instead of restOfDivision
+//NOTE: use kNF or kNF2 instead of restOfDivision
 inline ideal gcone::ffG(const ideal &H, const ideal &G)
 {
@@ -2356 +2304 @@
         for (int ii=0;ii<size;ii++)
         {
-                poly temp1;//=pInit();
-                poly temp2;//=pInit();
-                poly temp3;=pInit();//polys to temporarily store values for pSub
-                res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
-                temp1=pCopy(H->m[ii]);//TRY
-                temp2=pCopy(res->m[ii]);
-                NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
-                temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));//TRY
-                temp3=pSub(temp1, temp2);//TRY
-                temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0)));NOTRY
+//              poly temp1;//=pInit();
+//              poly temp2;//=pInit();
+                poly temp3;//=pInit();//polys to temporarily store values for pSub
+//              res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
+//              temp1=pCopy(H->m[ii]);//TRY
+//              temp2=pCopy(res->m[ii]);
+                //NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
+//              temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));//TRY
+//              temp3=pSub(temp1, temp2);//TRY
+                temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0)));//NOTRY
                 res->m[ii]=pCopy(temp3);
-                res->m[ii]=pSub(temp1,temp2); //buggy
-                cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
-                pDelete(&temp1);//TRY
-                pDelete(&temp2);
+                //res->m[ii]=pSub(temp1,temp2); //buggy         
+                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
+//              pDelete(&temp1);//TRY
+//              pDelete(&temp2);
                 pDelete(&temp3);
         }
         return res;
 }
-
+        
 /** \brief Preprocessing of inequlities
 * Do some preprocessing on the matrix of inequalities
@@ -2388 +2336 @@
         int num_elts=0;
         int offset=0;*/
-        Remove zeroes (and strictly pos rows?)
+        //Remove zeroes (and strictly pos rows?)
         for(int ii=0;ii<ddineq->rowsize;ii++)
         {
                 int64vec *iv = new int64vec(this->numVars);
-                int64vec *ivNull = new int64vec(this->numVars);Needed for intvec64::compare(*int64vec)
+                int64vec *ivNull = new int64vec(this->numVars);//Needed for intvec64::compare(*int64vec)
                 int posCtr=0;
                 for(int jj=0;jj<this->numVars;jj++)
                 {
                         (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
-                        if((*iv)[jj]>0)check for strictly pos rows
+                        if((*iv)[jj]>0)//check for strictly pos rows
                                 posCtr++;
-                        Behold! This will delete the row for the standard simplex!
+                        //Behold! This will delete the row for the standard simplex!
                 }
-                if( (iv->compare(0)==0) || (posCtr==iv->length()) )
-                if( (posCtr==iv->length()) || (iv->compare(ivNull)==0) )
+//              if( (iv->compare(0)==0) || (posCtr==iv->length()) )
+                if( (posCtr==iv->length()) || (iv->compare(ivNull)==0) )                
                 {
                         dd_MatrixRowRemove(&ddineq,ii+1);
-                        ii--;Yes. This is on purpose
+                        ii--;//Yes. This is on purpose
                 }
                 delete iv;
                 delete ivNull;
         }
-        Remove duplicates of rows
-        posRowsArray=NULL;
-        num_alloc=0;
-        num_elts=0;
-        offset=0;
-        int num_newRows = ddineq->rowsize;
-        for(int ii=0;ii<ddineq->rowsize-1;ii++)
-        for(int ii=0;ii<num_newRows-1;ii++)
-        {
-                int64vec *iv = new int64vec(this->numVars);//1st vector to check against
-                for(int jj=0;jj<this->numVars;jj++)
-                        (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
-                for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++)
-                {
-                        int64vec *ivCheck = new int64vec(this->numVars);//Checked against iv
-                        for(int kk=0;kk<this->numVars;kk++)
-                                (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
-                        if (iv->compare(ivCheck)==0)
-                        {
- //                             cout << "=" << endl;
- //                             num_alloc++;
- //                             void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int)));
- //                             if(!tmp)
- //                             {
- //                                     WerrorS("Woah dude! Couldn't realloc memory\n");
- //                                     exit(-1);
- //                             }
- //                             posRowsArray = (int*)tmp;
- //                             posRowsArray[num_elts]=jj;
- //                             num_elts++;
-                                dd_MatrixRowRemove(&ddineq,jj+1);
-                                num_newRows = ddineq->rowsize;
-                        }
-                        delete ivCheck;
-                }
-                delete iv;
-        }
-        for(int ii=0;ii<num_elts;ii++)
-        {
-                dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
-                offset++;
-        }
-        free(posRowsArray);
-        Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/
-}preprocessInequalities
-
+        //Remove duplicates of rows
+//      posRowsArray=NULL;
+//      num_alloc=0;
+//      num_elts=0;
+//      offset=0;
+//      int num_newRows = ddineq->rowsize;
+//      for(int ii=0;ii<ddineq->rowsize-1;ii++)
+//      for(int ii=0;ii<num_newRows-1;ii++)
+//      {
+//              int64vec *iv = new int64vec(this->numVars);//1st vector to check against
+//              for(int jj=0;jj<this->numVars;jj++)
+//                      (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
+//              for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++)
+//              {
+//                      int64vec *ivCheck = new int64vec(this->numVars);//Checked against iv
+//                      for(int kk=0;kk<this->numVars;kk++)
+//                              (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
+//                      if (iv->compare(ivCheck)==0)
+//                      {
+// //                           cout << "=" << endl;
+// //                           num_alloc++;
+// //                           void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int)));
+// //                           if(!tmp)
+// //                           {
+// //                                   WerrorS("Woah dude! Couldn't realloc memory\n");
+// //                                   exit(-1);
+// //                           }
+// //                           posRowsArray = (int*)tmp;
+// //                           posRowsArray[num_elts]=jj;
+// //                           num_elts++;
+//                              dd_MatrixRowRemove(&ddineq,jj+1);
+//                              num_newRows = ddineq->rowsize;
+//                      }
+//                      delete ivCheck;
+//              }
+//              delete iv;
+//      }
+//      for(int ii=0;ii<num_elts;ii++)
+//      {
+//              dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
+//              offset++;
+//      }
+//      free(posRowsArray);
+        //Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/        
+}//preprocessInequalities
+        
 /** \brief Compute a Groebner Basis
  *
@@ -2462 +2410 @@
  *\return void
  */
-inline void gcone::getGB(const ideal &inputIdeal)
+inline void gcone::getGB(const ideal &inputIdeal)               
 {
         BITSET save=test;
@@ -2471 +2419 @@
         idNorm(gb);
         idSkipZeroes(gb);
-        this->gcBasis=gb;       write the GB into gcBasis
+        this->gcBasis=gb;       //write the GB into gcBasis
         test=save;
-}void getGB
-
+}//void getGB
+                
 /** \brief Compute the negative of a given int64vec
-        */
+        */              
 static int64vec* ivNeg(/*const*/int64vec *iv)
-{       Hm, switching to int64vec const int64vec does no longer work
-        int64vec *res; = new int64vec(iv->length());
+{       //Hm, switching to int64vec const int64vec does no longer work
+        int64vec *res;// = new int64vec(iv->length());
         res=iv64Copy(iv);
-        *res *= (int)-1;
+        *res *= (int)-1;                                                
         return res;
 }
@@ -2489 +2437 @@
 *
 */
-static int dotProduct(const int64vec &iva, const int64vec &ivb)
-{
-        int res=0;
+static int dotProduct(const int64vec &iva, const int64vec &ivb)                         
+{                       
+        int res=0;      
         for (int i=0;i<pVariables;i++)
         {
- #ifndef NDEBUG
-        (const_cast<int64vec*>(&iva))->show(1,0); (const_cast<int64vec*>(&ivb))->show(1,0);
- #endif
+// #ifndef NDEBUG
+//      (const_cast<int64vec*>(&iva))->show(1,0); (const_cast<int64vec*>(&ivb))->show(1,0);
+// #endif
                 res = res+(iva[i]*ivb[i]);
         }
@@ -2506 +2454 @@
  */
 static bool isParallel(const int64vec &a,const int64vec &b)
-{
-/*#ifdef gfanp
+{       
+/*#ifdef gfanp  
         timeval start, end;
         gettimeofday(&start, 0);
-#endif*/
+#endif*/                
         bool res;
         int lhs=dotProduct(a,b)*dotProduct(a,b);
         int rhs=dotProduct(a,a)*dotProduct(b,b);
- #ifdef gfanp
-        gettimeofday(&end, 0);
-        t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
- #endif
-        return res;
+// #ifdef gfanp
+//      gettimeofday(&end, 0);
+//      t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
+// #endif       
+//      return res;
         return res = (lhs==rhs)?TRUE:FALSE;
-}bool isParallel
+}//bool isParallel
 
 /** \brief Compute an interior point of a given cone
- * Result will be written into int64vec iv.
+ * Result will be written into int64vec iv. 
  * Any rational point is automatically converted into an integer.
  */
-void gcone::interiorPoint( dd_MatrixPtr &M, int64vec &iv) no const &M here since we want to remove redundant rows
+void gcone::interiorPoint( dd_MatrixPtr &M, int64vec &iv) //no const &M here since we want to remove redundant rows
 {
         dd_LPPtr lp,lpInt;
@@ -2532 +2480 @@
         dd_LPSolverType solver=dd_DualSimplex;
         dd_LPSolutionPtr lpSol=NULL;
-        dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
-        dd_rowindex ddnewpos;
-        dd_NumberType numb;
-        M->representation=dd_Inequality;
-
-        NOTE: Make this n-dimensional!
-        dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
-
+//      dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
+//      dd_rowindex ddnewpos;
+        dd_NumberType numb;     
+        //M->representation=dd_Inequality;
+                        
+        //NOTE: Make this n-dimensional!
+        //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
+                                                                        
         /*NOTE: Leave the following line commented out!
         * Otherwise it will slow down computations a great deal
         * */
-        dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
-        if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
+//      dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
+        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
         dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
         int jj=1;
@@ -2550 +2498 @@
         {
                 dd_set_si(posRestr->matrix[ii][jj],1);
-                jj++;
+                jj++;                                                   
         }
         dd_MatrixAppendTo(&M,posRestr);
@@ -2558 +2506 @@
         if (lp==NULL){WerrorS("LP is NULL");}
 #ifdef gfan_DEBUG
-                        dd_WriteLP(stdout,lp);
-#endif
-
+//                      dd_WriteLP(stdout,lp);
+#endif  
+                                        
         lpInt=dd_MakeLPforInteriorFinding(lp);
         if (err!=dd_NoError){WerrorS("Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint");}
 #ifdef gfan_DEBUG
-                        dd_WriteLP(stdout,lpInt);
-#endif
-        dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
+//                      dd_WriteLP(stdout,lpInt);
+#endif                  
+//      dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
         if (err!=dd_NoError)
         {
                 WerrorS("Error during dd_FindRelativeInterior in gcone::interiorPoint");
                 dd_WriteErrorMessages(stdout, err);
-        }
-        dd_LPSolve(lpInt,solver,&err);  This will not result in a point from the relative interior
-        if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}
+        }                       
+        dd_LPSolve(lpInt,solver,&err);  //This will not result in a point from the relative interior
+//      if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}                                                                       
         lpSol=dd_CopyLPSolution(lpInt);
-        if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
+//      if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
 #ifdef gfan_DEBUG
         printf("Interior point: ");
@@ -2583 +2531 @@
         }
         printf("\n");
-#endif
-        NOTE The following strongly resembles parts of makeInt.
-        Maybe merge sometimes
+#endif                  
+        //NOTE The following strongly resembles parts of makeInt. 
+        //Maybe merge sometimes
         mpz_t kgV; mpz_init(kgV);
         mpz_set_str(kgV,"1",10);
@@ -2609 +2557 @@
         }
 #ifdef gfan_DEBUG
-                        iv.show();
-                        cout << endl;
+//                      iv.show();
+//                      cout << endl;
 #endif
         mpq_clear(qkgV);
         mpz_clear(tmp);
         mpz_clear(den);
-        mpz_clear(kgV);
-
+        mpz_clear(kgV);                 
+                        
         dd_FreeLPSolution(lpSol);
         dd_FreeLPData(lpInt);
         dd_FreeLPData(lp);
-        set_free(ddlinset);
-        set_free(ddredrows);
-
-}void interiorPoint(dd_MatrixPtr const &M)
+//      set_free(ddlinset);
+//      set_free(ddredrows);    
+                        
+}//void interiorPoint(dd_MatrixPtr const &M)
 
 /** Computes an interior point of a cone by taking two interior points a,b from two different facets
 * and then computing b+(a-b)/2
-* Of course this only works for flippable facets
-* Two cases may occur:
+* Of course this only works for flippable facets 
+* Two cases may occur: 
 * 1st: There are only two facets who share the only strictly positive ray
 * 2nd: There are at least two facets which have a distinct positive ray
@@ -2637 +2585 @@
 * positive => these lie in a plane and thus their sum is not from relative interior.
 * So let's just sum up all rays, find one strictly positive and shift the point along that ray
-*
+* 
 * Used by noRevS
 *NOTE no longer used nor maintained. MM Mar 9, 2010
 */
- void gcone::interiorPoint2()
- {//idPrint(this->gcBasis);
- #ifdef gfan_DEBUG
-        if(this->ivIntPt!=NULL)
-                WarnS("Interior point already exists - ovrewriting!");
- #endif
-        facet *f1 = this->facetPtr;
-        facet *f2 = NULL;
-        int64vec *intF1=NULL;
-        while(f1!=NULL)
-        {
-                if(f1->isFlippable)
-                {
-                        facet *f1Ray = f1->codim2Ptr;
-                        while(f1Ray!=NULL)
-                        {
-                                const int64vec *check=f1Ray->getRef2FacetNormal();
-                                if(iv64isStrictlyPositive(check))
-                                {
-                                        intF1=iv64Copy(check);
-                                        break;
-                                }
-                                f1Ray=f1Ray->next;
-                        }
-                }
-                if(intF1!=NULL)
-                        break;
-                f1=f1->next;
-        }
-        if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
-                f2=f1->next;
-        else
-                f2=this->facetPtr;
-        if(intF1==NULL && hasHomInput==TRUE)
-        {
-                intF1 = new int64vec(this->numVars);
-                for(int ii=0;ii<this->numVars;ii++)
-                        (*intF1)[ii]=1;
-        }
-        assert(f1); assert(f2);
-        int64vec *intF2=f2->getInteriorPoint();
-        mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
-        mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
-        mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
-        for(int ii=0;ii<this->numVars;ii++)
-        {
-                mpq_init(qPosRay[ii]);
-                mpq_init(qIntPt[ii]);
-                mpq_init(qPosIntPt[ii]);
-        }
-        //Compute a+((b-a)/2) && Convert intF1 to mpq
-        for(int ii=0;ii<this->numVars;ii++)
-        {
-                mpq_t a,b;
-                mpq_init(a); mpq_init(b);
-                mpq_set_si(a,(*intF1)[ii],1);
-                mpq_set_si(b,(*intF2)[ii],1);
-                mpq_t diff;
-                mpq_init(diff);
-                mpq_sub(diff,b,a);      //diff=b-a
-                mpq_t quot;
-                mpq_init(quot);
-                mpq_div_2exp(quot,diff,1);      //quot=diff/2=(b-a)/2
-                mpq_clear(diff);
-                //Don't be clever and reuse diff here
-                mpq_t sum; mpq_init(sum);
-                mpq_add(sum,b,quot);    //sum=b+quot=a+(b-a)/2
-                mpq_set(qIntPt[ii],sum);
-                mpq_clear(sum);
-                mpq_clear(quot);
-                mpq_clear(a); mpq_clear(b);
-                //Now for intF1
-                mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
-        }
-        //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
-        while(TRUE)
-        {
-                bool success=FALSE;
-                int posCtr=0;
-                for(int ii=0;ii<this->numVars;ii++)
-                {
-                        mpq_t sum; mpq_init(sum);
-                        mpq_add(sum,qPosRay[ii],qIntPt[ii]);
-                        mpq_set(qPosIntPt[ii],sum);
-                        mpq_clear(sum);
-                        if(mpq_sgn(qPosIntPt[ii])==1)
-                                posCtr++;
-                }
-                if(posCtr==this->numVars)//qPosIntPt > 0
-                        break;
-                else
-                {
-                        mpq_t qTwo; mpq_init(qTwo);
-                        mpq_set_ui(qTwo,2,1);
-                        for(int jj=0;jj<this->numVars;jj++)
-                        {
-                                mpq_t tmp; mpq_init(tmp);
-                                mpq_mul(tmp,qPosRay[jj],qTwo);
-                                mpq_set( qPosRay[jj], tmp);
-                                mpq_clear(tmp);
-                        }
-                        mpq_clear(qTwo);
-                }
-        }//while
-        //Now qPosIntPt ought to be >0, so convert back to int :D
-        /*Compute lcm of the denominators*/
-        mpz_t *denom = new mpz_t[this->numVars];
-        mpz_t tmp,kgV;
-        mpz_init(tmp); mpz_init(kgV);
-        for (int ii=0;ii<this->numVars;ii++)
-        {
-                mpz_t z;
-                mpz_init(z);
-                mpq_get_den(z,qPosIntPt[ii]);
-                mpz_init(denom[ii]);
-                mpz_set( denom[ii], z);
-                mpz_clear(z);
-        }
-
-        mpz_set(tmp,denom[0]);
-        for (int ii=0;ii<this->numVars;ii++)
-        {
-                mpz_lcm(kgV,tmp,denom[ii]);
-                mpz_set(tmp,kgV);
-        }
-        mpz_clear(tmp);
-        /*Multiply the nominators by kgV*/
-        mpq_t qkgV,res;
-        mpq_init(qkgV);
-        mpq_canonicalize(qkgV);
-        mpq_init(res);
-        mpq_canonicalize(res);
-
-        mpq_set_num(qkgV,kgV);
-        int64vec *n=new int64vec(this->numVars);
-        for (int ii=0;ii<this->numVars;ii++)
-        {
-                mpq_canonicalize(qPosIntPt[ii]);
-                mpq_mul(res,qkgV,qPosIntPt[ii]);
-                (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
-        }
-        this->setIntPoint(n);
-        delete n;
-        delete [] qPosIntPt;
-        delete [] denom;
-        delete [] qPosRay;
-        delete [] qIntPt;
-        mpz_clear(kgV);
-        mpq_clear(qkgV); mpq_clear(res);
- }
-
+// void gcone::interiorPoint2()
+// {//idPrint(this->gcBasis);
+// #ifdef gfan_DEBUG
+//      if(this->ivIntPt!=NULL)
+//              WarnS("Interior point already exists - ovrewriting!");
+// #endif
+//      facet *f1 = this->facetPtr;
+//      facet *f2 = NULL;
+//      int64vec *intF1=NULL;
+//      while(f1!=NULL)
+//      {
+//              if(f1->isFlippable)
+//              {
+//                      facet *f1Ray = f1->codim2Ptr;
+//                      while(f1Ray!=NULL)
+//                      {
+//                              const int64vec *check=f1Ray->getRef2FacetNormal();
+//                              if(iv64isStrictlyPositive(check))
+//                              {
+//                                      intF1=iv64Copy(check);
+//                                      break;
+//                              }                               
+//                              f1Ray=f1Ray->next;
+//                      }
+//              }
+//              if(intF1!=NULL)
+//                      break;
+//              f1=f1->next;
+//      }
+//      if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
+//              f2=f1->next;
+//      else
+//              f2=this->facetPtr;
+//      if(intF1==NULL && hasHomInput==TRUE)
+//      {
+//              intF1 = new int64vec(this->numVars);
+//              for(int ii=0;ii<this->numVars;ii++)
+//                      (*intF1)[ii]=1;
+//      }
+//      assert(f1); assert(f2);
+//      int64vec *intF2=f2->getInteriorPoint();
+//      mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
+//      mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
+//      mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
+//      for(int ii=0;ii<this->numVars;ii++)
+//      {
+//              mpq_init(qPosRay[ii]);
+//              mpq_init(qIntPt[ii]);
+//              mpq_init(qPosIntPt[ii]);
+//      }       
+//      //Compute a+((b-a)/2) && Convert intF1 to mpq
+//      for(int ii=0;ii<this->numVars;ii++)
+//      {
+//              mpq_t a,b;
+//              mpq_init(a); mpq_init(b);
+//              mpq_set_si(a,(*intF1)[ii],1);
+//              mpq_set_si(b,(*intF2)[ii],1);
+//              mpq_t diff;
+//              mpq_init(diff);
+//              mpq_sub(diff,b,a);      //diff=b-a
+//              mpq_t quot;
+//              mpq_init(quot);
+//              mpq_div_2exp(quot,diff,1);      //quot=diff/2=(b-a)/2
+//              mpq_clear(diff);
+//              //Don't be clever and reuse diff here
+//              mpq_t sum; mpq_init(sum);
+//              mpq_add(sum,b,quot);    //sum=b+quot=a+(b-a)/2
+//              mpq_set(qIntPt[ii],sum);
+//              mpq_clear(sum);
+//              mpq_clear(quot);
+//              mpq_clear(a); mpq_clear(b);
+//              //Now for intF1
+//              mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
+//      }
+//      //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
+//      while(TRUE)
+//      {       
+//              bool success=FALSE;
+//              int posCtr=0;   
+//              for(int ii=0;ii<this->numVars;ii++)
+//              {
+//                      mpq_t sum; mpq_init(sum);
+//                      mpq_add(sum,qPosRay[ii],qIntPt[ii]);
+//                      mpq_set(qPosIntPt[ii],sum);
+//                      mpq_clear(sum);
+//                      if(mpq_sgn(qPosIntPt[ii])==1)
+//                              posCtr++;
+//              }
+//              if(posCtr==this->numVars)//qPosIntPt > 0 
+//                      break;
+//              else
+//              {
+//                      mpq_t qTwo; mpq_init(qTwo);
+//                      mpq_set_ui(qTwo,2,1);
+//                      for(int jj=0;jj<this->numVars;jj++)
+//                      {
+//                              mpq_t tmp; mpq_init(tmp);
+//                              mpq_mul(tmp,qPosRay[jj],qTwo);
+//                              mpq_set( qPosRay[jj], tmp);
+//                              mpq_clear(tmp);
+//                      }
+//                      mpq_clear(qTwo);
+//              }
+//      }//while
+//      //Now qPosIntPt ought to be >0, so convert back to int :D
+//      /*Compute lcm of the denominators*/
+//      mpz_t *denom = new mpz_t[this->numVars];
+//      mpz_t tmp,kgV;
+//      mpz_init(tmp); mpz_init(kgV);
+//      for (int ii=0;ii<this->numVars;ii++)
+//      {
+//              mpz_t z;
+//              mpz_init(z);
+//              mpq_get_den(z,qPosIntPt[ii]);
+//              mpz_init(denom[ii]);
+//              mpz_set( denom[ii], z);
+//              mpz_clear(z);                           
+//      }
+//              
+//      mpz_set(tmp,denom[0]);
+//      for (int ii=0;ii<this->numVars;ii++)
+//      {
+//              mpz_lcm(kgV,tmp,denom[ii]);
+//              mpz_set(tmp,kgV);                               
+//      }
+//      mpz_clear(tmp); 
+//      /*Multiply the nominators by kgV*/
+//      mpq_t qkgV,res;
+//      mpq_init(qkgV);
+//      mpq_canonicalize(qkgV);         
+//      mpq_init(res);
+//      mpq_canonicalize(res);
+//                              
+//      mpq_set_num(qkgV,kgV);
+//      int64vec *n=new int64vec(this->numVars);
+//      for (int ii=0;ii<this->numVars;ii++)
+//      {
+//              mpq_canonicalize(qPosIntPt[ii]);
+//              mpq_mul(res,qkgV,qPosIntPt[ii]);
+//              (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
+//      }
+//      this->setIntPoint(n);
+//      delete n;
+//      delete [] qPosIntPt;
+//      delete [] denom;
+//      delete [] qPosRay;
+//      delete [] qIntPt;
+//      mpz_clear(kgV);
+//      mpq_clear(qkgV); mpq_clear(res);
+// }
+        
 /** \brief Copy a ring and add a weighted ordering in first place
- *
+ * 
  */
-ring gcone::rCopyAndAddWeight(const ring &r, int64vec *ivw)
+ring gcone::rCopyAndAddWeight(const ring &r, int64vec *ivw)                             
 {
         ring res=rCopy0(r);
         int jj;
-
+                        
         omFree(res->order);
         res->order =(int *)omAlloc0(4*sizeof(int/*64*/));
@@ -2808 +2756 @@
         omfree(res->wvhdl);
         res->wvhdl =(int **)omAlloc0(4*sizeof(int/*64*/**));
-
+                        
         res->order[0]=ringorder_a/*64*/;
         res->block0[0]=1;
         res->block1[0]=res->N;
-        res->order[1]=ringorder_dp;     basically useless, since that should never be used
+        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                    
         res->block0[1]=1;
         res->block1[1]=res->N;
@@ -2820 +2768 @@
         int/*64*/ *A=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
         for (jj=0;jj<length;jj++)
-        {
+        {                               
                 A[jj]=(*ivw)[jj];
                 if((*ivw)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::rCopyAndAddWeight!\n");
-        }
+        }                       
         res->wvhdl[0]=(int *)A;
         res->block1[0]=length;
-
+                        
         rComplete(res);
         return res;
-}rCopyAndAdd
-
-ring gcone::rCopyAndAddWeight2(const ring &r,const int64vec *ivw, const int64vec *fNormal)
+}//rCopyAndAdd
+                
+ring gcone::rCopyAndAddWeight2(const ring &r,const int64vec *ivw, const int64vec *fNormal)                              
 {
         ring res=rCopy0(r);
-
+                        
         omFree(res->order);
         res->order =(int *)omAlloc0(5*sizeof(int/*64*/));
@@ -2843 +2791 @@
         omfree(res->wvhdl);
         res->wvhdl =(int **)omAlloc0(5*sizeof(int/*64*/**));
-
+                        
         res->order[0]=ringorder_a/*64*/;
         res->block0[0]=1;
@@ -2850 +2798 @@
         res->block0[1]=1;
         res->block1[1]=res->N;
-
+        
         res->order[2]=ringorder_dp;
         res->block0[2]=1;
         res->block1[2]=res->N;
-
+        
         res->order[3]=ringorder_C;
 
@@ -2861 +2809 @@
         int/*64*/ *A2=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
         for (int jj=0;jj<length;jj++)
-        {
+        {                               
                 A1[jj]=(*ivw)[jj];
                 A2[jj]=-(*fNormal)[jj];
                 if((*ivw)[jj]>=INT_MAX || (*fNormal)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::rCopyAndAddWeight2!\n");
-        }
+        }                       
         res->wvhdl[0]=(int *)A1;
         res->block1[0]=length;
@@ -2873 +2821 @@
         return res;
 }
-
-NOTE not needed anywhere
- ring rCopyAndChangeWeight(ring const &r, int64vec *ivw)
- {
-        ring res=rCopy0(currRing);
-        rComplete(res);
-        rSetWeightVec(res,(int64*)ivw);
-        //rChangeCurrRing(rnew);
-        return res;
- }
-
+                
+//NOTE not needed anywhere
+// ring rCopyAndChangeWeight(ring const &r, int64vec *ivw)
+// {
+//      ring res=rCopy0(currRing);
+//      rComplete(res);
+//      rSetWeightVec(res,(int64*)ivw);
+//      //rChangeCurrRing(rnew);
+//      return res;
+// }
+                
 /** \brief Checks whether a given facet is a search facet
  * Determines whether a given facet of a cone is the search facet of a neighbouring cone
@@ -2890 +2838 @@
  * 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.
+ * If this is the case, then our facet is indeed a search facet and TRUE is retuned. 
 */
-removed with r12286
-
+//removed with r12286
+                
 /** \brief Check for equality of two intvecs
  */
@@ -2909 +2857 @@
         return res;
 }
-
+                
 /** \brief The reverse search algorithm
  */
-removed with r12286
+//removed with r12286
 /** \brief Compute the lineality/homogeneity space
 * It is the kernel of the inequality matrix Ax=0
@@ -2921 +2869 @@
         dd_MatrixPtr res;
         dd_MatrixPtr ddineq;
-        ddineq=dd_CopyMatrix(this->ddFacets);
-        Add a row of 0s in 0th place
+        ddineq=dd_CopyMatrix(this->ddFacets);   
+        //Add a row of 0s in 0th place
         dd_MatrixPtr ddAppendRowOfZeroes=dd_CreateMatrix(1,this->numVars+1);
         dd_MatrixPtr ddFoo=dd_AppendMatrix(ddAppendRowOfZeroes,ddineq);
@@ -2929 +2877 @@
         ddineq=dd_CopyMatrix(ddFoo);
         dd_FreeMatrix(ddFoo);
-        Cohen starts here
-        int dimKer=0;Cohen calls this r
-        int m=ddineq->rowsize;Rows
-        int n=ddineq->colsize;Cols
+        //Cohen starts here
+        int dimKer=0;//Cohen calls this r
+        int m=ddineq->rowsize;//Rows
+        int n=ddineq->colsize;//Cols
         int c[m+1];
         int d[n+1];
@@ -2938 +2886 @@
                 c[ii]=0;
         for(int ii=0;ii<n;ii++)
-                d[ii]=0;
-
+                d[ii]=0;        
+        
         for(int k=1;k<n;k++)
         {
-                //Let's find a j s.t. m[j][k]!=0 && c[j]=0
-                int condCtr=0;Check each row for zeroness
+                //Let's find a j s.t. m[j][k]!=0 && c[j]=0              
+                int condCtr=0;//Check each row for zeroness
                 for(int j=1;j<m;j++)
                 {
@@ -2955 +2903 @@
                                 mpq_set(ratd,quot);
                                 mpq_canonicalize(ratd);
-
+                
                                 mpq_set_str(ddineq->matrix[j][k],"-1",10);
                                 for(int ss=k+1;ss<n;ss++)
                                 {
                                         mpq_t prod; mpq_init(prod);
-                                        mpq_mul(prod, ratd, ddineq->matrix[j][ss]);
+                                        mpq_mul(prod, ratd, ddineq->matrix[j][ss]);                             
                                         mpq_set(ddineq->matrix[j][ss],prod);
                                         mpq_canonicalize(ddineq->matrix[j][ss]);
                                         mpq_clear(prod);
-                                }
+                                }               
                                 for(int ii=1;ii<m;ii++)
                                 {
@@ -2970 +2918 @@
                                         {
                                                 mpq_set(ratd,ddineq->matrix[ii][k]);
-                                                mpq_set_str(ddineq->matrix[ii][k],"0",10);
+                                                mpq_set_str(ddineq->matrix[ii][k],"0",10);                      
                                                 for(int ss=k+1;ss<n;ss++)
                                                 {
@@ -2984 +2932 @@
                                         }
                                 }
-                                c[j]=k;
+                                c[j]=k;         
                                 d[k]=j;
-                                mpq_clear(quot); mpq_clear(ratd); mpq_clear(one);
+                                mpq_clear(quot); mpq_clear(ratd); mpq_clear(one);       
                         }
                         else
                                 condCtr++;
-                }
-                if(condCtr==m-1)        Why -1 ???
+                }                       
+                if(condCtr==m-1)        //Why -1 ???
                 {
                         dimKer++;
                         d[k]=0;
-                        break;//goto _4;
-                }else{
+//                      break;//goto _4;
+                }//else{
                 /*Eliminate*/
-        }for(k
-        /*Output kernel, i.e. set gcone::dd_LinealitySpace here*/
-        gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
+        }//for(k
+        /*Output kernel, i.e. set gcone::dd_LinealitySpace here*/       
+//      gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
         res = dd_CreateMatrix(dimKer,this->numVars+1);
         int row=-1;
@@ -3012 +2960 @@
                                 if(d[ii]>0)
                                         mpq_set(res->matrix[row][ii],ddineq->matrix[d[ii]][kk]);
-                                else if(ii==kk)
+                                else if(ii==kk)                         
                                         mpq_set_str(res->matrix[row][ii],"1",10);
                                 mpq_canonicalize(res->matrix[row][ii]);
@@ -3020 +2968 @@
         dd_FreeMatrix(ddineq);
         return(res);
-        Better safe than sorry:
-        NOTE dd_LinealitySpace contains RATIONAL ENTRIES
-        Therefore if you need integer ones: CALL gcone::makeInt(...) method
-}
+        //Better safe than sorry:
+        //NOTE dd_LinealitySpace contains RATIONAL ENTRIES
+        //Therefore if you need integer ones: CALL gcone::makeInt(...) method
+}       
 
 
@@ -3034 +2982 @@
  */
 void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
-{
-        facet *SearchListRoot = new facet(); The list containing ALL facets we come across
+{       
+        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
         facet *SearchListAct;
         SearchListAct = NULL;
-        SearchListAct = SearchListRoot;
+        //SearchListAct = SearchListRoot;                       
         gcone *gcAct;
         gcAct = &gcRoot;
-        gcone *gcPtr;   Pointer to end of linked list of cones
+        gcone *gcPtr;   //Pointer to end of linked list of cones
         gcPtr = &gcRoot;
-        gcone *gcNext;  Pointer to next cone to be visited
+        gcone *gcNext;  //Pointer to next cone to be visited
         gcNext = &gcRoot;
         gcone *gcHead;
-        gcHead = &gcRoot;
+        gcHead = &gcRoot;                       
         facet *fAct;
-        fAct = gcAct->facetPtr;
+        fAct = gcAct->facetPtr;                         
         ring rAct;
         rAct = currRing;
-        int UCNcounter=gcAct->getUCN();
+        int UCNcounter=gcAct->getUCN(); 
 #ifdef gfan_DEBUG
         printf("NoRevs\n");
         printf("Facets are:\n");
         gcAct->showFacets();
-#endif
+#endif                  
         /*Compute lineality space here
         * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/
         if(dd_LinealitySpace==NULL)
                 dd_LinealitySpace = gcAct->computeLinealitySpace();
-        /*End of lineality space computation*/
-        gcAct->getCodim2Normals(*gcAct);
+        /*End of lineality space computation*/          
+        //gcAct->getCodim2Normals(*gcAct);
         if(fAct->codim2Ptr==NULL)
-                gcAct->getExtremalRays(*gcAct);
+                gcAct->getExtremalRays(*gcAct);         
         /* Make a copy of the facet list of first cone
            Since the operations getCodim2Normals and normalize affect the facets
-           we must not memcpy them before these ops! */
-        /*facet *codim2Act; codim2Act = NULL;
+           we must not memcpy them before these ops! */ 
+        /*facet *codim2Act; codim2Act = NULL;                   
         facet *sl2Root;
-        facet *sl2Act;*/
+        facet *sl2Act;*/                        
         for(int ii=0;ii<this->numFacets;ii++)
         {
-                only copy flippable facets! NOTE: We do not need flipRing or any such information.
+                //only copy flippable facets! NOTE: We do not need flipRing or any such information.
                 if(fAct->isFlippable==TRUE)
                 {
                         /*Using shallow copy here*/
 #ifdef SHALLOW
-                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) 1st facet may be non-flippable
+                        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!
+                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
                         }
                         else
-                        {
+                        {                       
                                 SearchListAct->next = fAct->shallowCopy(*fAct);
-                                SearchListAct = SearchListAct->next;
+                                SearchListAct = SearchListAct->next;                                            
                         }
                         fAct=fAct->next;
@@ -3094 +3042 @@
                         int64vec *fNormal;
                         fNormal = fAct->getFacetNormal();
-                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) 1st facet may be non-flippable
-                        {
-                                SearchListAct = SearchListRoot; SearchListRoot is already 'new'ed at beginning of method!
+                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
+                        {                                               
+                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
                         }
                         else
-                        {
+                        {                       
                                 SearchListAct->next = new facet();
-                                SearchListAct = SearchListAct->next;
+                                SearchListAct = SearchListAct->next;                                            
                         }
                         SearchListAct->setFacetNormal(fNormal);
@@ -3107 +3055 @@
                         SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
                         SearchListAct->isFlippable=TRUE;
-                        Copy int point as well
+                        //Copy int point as well
                         int64vec *ivIntPt;
                         ivIntPt = fAct->getInteriorPoint();
                         SearchListAct->setInteriorPoint(ivIntPt);
                         delete(ivIntPt);
-                        Copy codim2-facets
+                        //Copy codim2-facets
                         facet *codim2Act;
-                        codim2Act = NULL;
+                        codim2Act = NULL;                       
                         facet *sl2Root;
-                        facet *sl2Act;
+                        facet *sl2Act;                  
                         codim2Act=fAct->codim2Ptr;
                         SearchListAct->codim2Ptr = new facet(2);
                         sl2Root = SearchListAct->codim2Ptr;
-                        sl2Act = sl2Root;
+                        sl2Act = sl2Root;                       
                         for(int jj=0;jj<fAct->numCodim2Facets;jj++)
-                        for(int jj=0;jj<fAct->numRays-1;jj++)
+//                      for(int jj=0;jj<fAct->numRays-1;jj++)
                         {
                                 int64vec *f2Normal;
                                 f2Normal = codim2Act->getFacetNormal();
                                 if(jj==0)
-                                {
+                                {                                               
                                         sl2Act = sl2Root;
                                         sl2Act->setFacetNormal(f2Normal);
@@ -3137 +3085 @@
                                         sl2Act->setFacetNormal(f2Normal);
                                 }
-                                delete f2Normal;
+                                delete f2Normal;                                
                                 codim2Act = codim2Act->next;
                         }
@@ -3143 +3091 @@
                         delete fNormal;
 #endif
-                }if(fAct->isFlippable==TRUE)
+                }//if(fAct->isFlippable==TRUE)                  
                 else {fAct = fAct->next;}
-        }End of copying facets into SLA
-
-        SearchListAct = SearchListRoot; Set to beginning of list
+        }//End of copying facets into SLA
+        
+        SearchListAct = SearchListRoot; //Set to beginning of list
         /*Make SearchList doubly linked*/
 #ifndef NDEBUG
@@ -3164 +3112 @@
                 if(SearchListAct->next!=NULL){
 #endif
-                        SearchListAct->next->prev = SearchListAct;
+                        SearchListAct->next->prev = SearchListAct;                                      
                 }
                 SearchListAct = SearchListAct->next;
         }
-        SearchListAct = SearchListRoot; Set to beginning of List
-
-        fAct = gcAct->facetPtr;//???
-        gcAct->writeConeToFile(*gcAct);
+        SearchListAct = SearchListRoot; //Set to beginning of List
+        
+//      fAct = gcAct->facetPtr;//???
+        gcAct->writeConeToFile(*gcAct);                 
         /*End of initialisation*/
-
+        
         fAct = SearchListAct;
         /*2nd step
           Choose a facet from SearchList, flip it and forget the previous cone
           We always choose the first facet from SearchList as facet to be flipped
-        */
-        while( (SearchListAct!=NULL))&& counter<490)
-        {NOTE See to it that the cone is only changed after ALL facets have been flipped!
-                fAct = SearchListAct;
+        */      
+        while( (SearchListAct!=NULL))//&& counter<490)
+        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
+                fAct = SearchListAct;           
                 while(fAct!=NULL)
-                while( (fAct->getUCN() == fAct->next->getUCN()) )
-                {       Since SLA should only contain flippables there should be no need to check for that
+//              while( (fAct->getUCN() == fAct->next->getUCN()) )               
+                {       //Since SLA should only contain flippables there should be no need to check for that                    
                         gcAct->flip2(gcAct->gcBasis,fAct);
-                        NOTE rCopy needed?
+                        //NOTE rCopy needed?
                         ring rTmp=rCopy(fAct->flipRing);
                         rComplete(rTmp);
                         rChangeCurrRing(rTmp);
-                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);copy constructor!
+                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);//copy constructor!
                         /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
                          * Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
@@ -3200 +3148 @@
                         */
                         idDelete((ideal *)&fAct->flipGB);
-                        rDelete(fAct->flipRing);
+                        rDelete(fAct->flipRing);                        
                         gcTmp->getConeNormals(gcTmp->gcBasis/*, FALSE*/);
-                        gcTmp->getCodim2Normals(*gcTmp);
+//                      gcTmp->getCodim2Normals(*gcTmp);
                         gcTmp->getExtremalRays(*gcTmp);
-                        NOTE Order rays lex here
-                        gcTmp->orderRays();
-                        //NOTE If flip2 is used we need to get an interior point of gcTmp
-                        // and replace gcTmp->baseRing with an appropriate ring with only
-                        // one weight
-                        gcTmp->interiorPoint2();
+                        //NOTE Order rays lex here
+                        gcTmp->orderRays();                     
+//                      //NOTE If flip2 is used we need to get an interior point of gcTmp
+//                      // and replace gcTmp->baseRing with an appropriate ring with only
+//                      // one weight
+//                      gcTmp->interiorPoint2();
                         gcTmp->replaceDouble_ringorder_a_ByASingleOne();
-                        gcTmp->ddFacets should not be needed anymore, so
+                        //gcTmp->ddFacets should not be needed anymore, so
                         dd_FreeMatrix(gcTmp->ddFacets);
 #ifdef gfan_DEBUG
-                        gcTmp->showFacets(1);
+//                      gcTmp->showFacets(1);
 #endif
                         /*add facets to SLA here*/
@@ -3228 +3176 @@
   #endif
                         SearchListRoot=tmp;
-                        SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
-#else
+                        //SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
+#else 
                         SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
-#endif ifdef SHALLOW
-                        gcTmp->writeConeToFile(*gcTmp);
+#endif //ifdef SHALLOW
+//                      gcTmp->writeConeToFile(*gcTmp);
                         if(gfanHeuristic==1)
                         {
                                 gcTmp->writeConeToFile(*gcTmp);
-                                idDelete((ideal*)&gcTmp->gcBasis);Whonder why?
+                                idDelete((ideal*)&gcTmp->gcBasis);//Whonder why?
                                 rDelete(gcTmp->baseRing);
-                        }
+                        }                       
 #ifdef gfan_DEBUG
                         if(SearchListRoot!=NULL)
                                 showSLA(*SearchListRoot);
-#endif
+#endif                  
                         rChangeCurrRing(gcAct->baseRing);
                         rDelete(rTmp);
-                        doubly linked for easier removal
+                        //doubly linked for easier removal
                         gcTmp->prev = gcPtr;
                         gcPtr->next=gcTmp;
                         gcPtr=gcPtr->next;
-                        Cleverly disguised exit condition follows
+                        //Cleverly disguised exit condition follows
                         if(fAct->getUCN() == fAct->next->getUCN())
                         {
                                 printf("Switching UCN from %i to %i\n",fAct->getUCN(),fAct->next->getUCN());
-                                fAct=fAct->next;
+                                fAct=fAct->next;                                
                         }
                         else
                         {
-                                rDelete(gcAct->baseRing);
-                                printf("break\n");
+                                //rDelete(gcAct->baseRing);
+//                              printf("break\n");
                                 break;
                         }
-                        fAct=fAct->next;
-                }while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
-                Search for cone with smallest UCN
+//                      fAct=fAct->next;
+                }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );          
+                //Search for cone with smallest UCN
 #ifndef NDEBUG
-  #if SIZEOF_LONG==8    64 bit
+  #if SIZEOF_LONG==8    //64 bit
                 while(gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && SearchListRoot!=NULL)
   #elif SIZEOF_LONG == 4
@@ -3272 +3220 @@
 #endif
 #ifdef NDEBUG
-                while(gcNext!=NULL && SearchListRoot!=NULL)
-#endif
-                {
+                while(gcNext!=NULL && SearchListRoot!=NULL)     
+#endif
+                {                       
                         if( gcNext->getUCN() == SearchListRoot->getUCN() )
                         {
@@ -3280 +3228 @@
                                 {
                                         gcAct = gcNext;
-                                        Seems better to not to use rCopy here
-                                        rAct=rCopy(gcAct->baseRing);
+                                        //Seems better to not to use rCopy here
+//                                      rAct=rCopy(gcAct->baseRing);
                                         rAct=gcAct->baseRing;
                                         rComplete(rAct);
-                                        rChangeCurrRing(rAct);
+                                        rChangeCurrRing(rAct);                                          
                                         break;
                                 }
@@ -3290 +3238 @@
                                 {
                                         gcone *gcDel;
-                                        gcDel = gcAct;
+                                        gcDel = gcAct;                                  
                                         gcAct = gcNext;
-                                        Read st00f from file &
-                                        implant the GB into gcAct
+                                        //Read st00f from file &
+                                        //implant the GB into gcAct
                                         readConeFromFile(gcAct->getUCN(), gcAct);
-                                        Kill the baseRing but ONLY if it is not the ring the computation started in!
-                                        if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
-                                                rDelete(gcDel->baseRing);
-                                        rAct=rCopy(gcAct->baseRing);
-                                        /*The ring change occurs in readConeFromFile, so as to
+                                        //Kill the baseRing but ONLY if it is not the ring the computation started in!
+//                                      if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
+//                                              rDelete(gcDel->baseRing);
+//                                      rAct=rCopy(gcAct->baseRing);
+                                        /*The ring change occurs in readConeFromFile, so as to 
                                         assure that gcAct->gcBasis belongs to the right ring*/
                                         rAct=gcAct->baseRing;
@@ -3305 +3253 @@
                                         rChangeCurrRing(rAct);
                                         break;
-                                }
+                                }                               
                         }
                         else if(gcNext->getUCN() < SearchListRoot->getUCN() )
                         {
-                                idDelete( (ideal*)&gcNext->gcBasis );
-                                rDelete(gcNext->baseRing);//TODO Why does this crash?
+                                idDelete( (ideal*)&gcNext->gcBasis );                           
+//                              rDelete(gcNext->baseRing);//TODO Why does this crash?
                         }
                         /*else
@@ -3320 +3268 @@
                                         if(gcDel->getUCN()!=1)
                                                 rDelete(gcDel->baseRing);
-                                }
-                        }*/
+                                }                                       
+                        }*/                     
                         gcNext = gcNext->next;
                 }
                 UCNcounter++;
-                SearchListAct = SearchListRoot;
+                SearchListAct = SearchListRoot;         
         }
         printf("\nFound %i cones - terminating\n", counter);
-}void noRevS(gcone &gc)
-
-
+}//void noRevS(gcone &gc)       
+                
+                
 /** \brief Make a set of rational vectors into integers
  *
@@ -3339 +3287 @@
  */
 void gcone::makeInt(const dd_MatrixPtr &M, const int line, int64vec &n)
-{
+{                       
         mpz_t *denom = new mpz_t[this->numVars];
         for(int ii=0;ii<this->numVars;ii++)
@@ -3345 +3293 @@
                 mpz_init_set_str(denom[ii],"0",10);
         }
-
+                
         mpz_t kgV,tmp;
         mpz_init(kgV);
         mpz_init(tmp);
-
+                        
         for (int ii=0;ii<(M->colsize)-1;ii++)
         {
@@ -3356 +3304 @@
                 mpq_get_den(z,M->matrix[line-1][ii+1]);
                 mpz_set( denom[ii], z);
-                mpz_clear(z);
-        }
-
+                mpz_clear(z);                           
+        }
+                                                
         /*Compute lcm of the denominators*/
         mpz_set(tmp,denom[0]);
@@ -3364 +3312 @@
         {
                 mpz_lcm(kgV,tmp,denom[ii]);
-                mpz_set(tmp,kgV);
-        }
-        mpz_clear(tmp);
+                mpz_set(tmp,kgV);                               
+        }
+        mpz_clear(tmp); 
         /*Multiply the nominators by kgV*/
         mpq_t qkgV,res;
@@ -3372 +3320 @@
         mpq_set_str(qkgV,"1",10);
         mpq_canonicalize(qkgV);
-
+                        
         mpq_init(res);
         mpq_set_str(res,"1",10);
         mpq_canonicalize(res);
-
+                        
         mpq_set_num(qkgV,kgV);
-
-        mpq_canonicalize(qkgV);
-        int ggT=1;
+                        
+//      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]=(int64)mpz_get_d(mpq_numref(res));
-                ggT=int64gcd(ggT,n[ii]);
+//              ggT=int64gcd(ggT,n[ii]);
         }
         int64 ggT=n[0];
         for(int ii=0;ii<this->numVars;ii++)
-                ggT=int64gcd(ggT,n[ii]);
-        Normalisation
+                ggT=int64gcd(ggT,n[ii]);        
+        //Normalisation
         if(ggT>1)
         {
@@ -3398 +3346 @@
         delete [] denom;
         mpz_clear(kgV);
-        mpq_clear(qkgV); mpq_clear(res);
-
+        mpq_clear(qkgV); mpq_clear(res);                        
+                        
 }
 /**
- * For all codim-2-facets we compute the gcd of the components of the facet normal and
+ * For all codim-2-facets we compute the gcd of the components of the facet normal and 
  * divide it out. Thus we get a normalized representation of each
  * (codim-2)-facet normal, i.e. a primitive vector
  * Actually we now also normalize the facet normals.
  */
- void gcone::normalize()
- {
-        int *ggT = new int;
-                *ggT=1;
-        facet *fAct;
-        facet *codim2Act;
-        fAct = this->facetPtr;
-        codim2Act = fAct->codim2Ptr;
-        while(fAct!=NULL)
-        {
-                int64vec *fNormal;
-                fNormal = fAct->getFacetNormal();
-                int *ggT = new int;
-                *ggT=1;
-                for(int ii=0;ii<this->numVars;ii++)
-                {
-                        *ggT=intgcd((*ggT),(*fNormal)[ii]);
-                }
-                if(*ggT>1)//We only need to do this if the ggT is non-trivial
-                {
- //                     int64vec *fCopy = fAct->getFacetNormal();
-                        for(int ii=0;ii<this->numVars;ii++)
-                                (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT);
-                        fAct->setFacetNormal(fNormal);
-                }
-                delete fNormal;
-                delete ggT;
-                /*And now for the codim2*/
-                while(codim2Act!=NULL)
-                {
-                        int64vec *n;
-                        n=codim2Act->getFacetNormal();
-                        int *ggT=new int;
-                        *ggT=1;
-                        for(int ii=0;ii<this->numVars;ii++)
-                        {
-                                *ggT = intgcd((*ggT),(*n)[ii]);
-                        }
-                        if(*ggT>1)
-                        {
-                                for(int ii=0;ii<this->numVars;ii++)
-                                {
-                                        (*n)[ii] = ((*n)[ii])/(*ggT);
-                                }
-                                codim2Act->setFacetNormal(n);
-                        }
-                        codim2Act = codim2Act->next;
-                        delete n;
-                        delete ggT;
-                }
-                fAct = fAct->next;
-        }
- }
-
-/** \brief Enqueue new facets into the searchlist
+// void gcone::normalize()
+// {
+//      int *ggT = new int;
+//              *ggT=1;
+//      facet *fAct;
+//      facet *codim2Act;
+//      fAct = this->facetPtr;
+//      codim2Act = fAct->codim2Ptr;
+//      while(fAct!=NULL)
+//      {
+//              int64vec *fNormal;
+//              fNormal = fAct->getFacetNormal();
+//              int *ggT = new int;
+//              *ggT=1;
+//              for(int ii=0;ii<this->numVars;ii++)
+//              {
+//                      *ggT=intgcd((*ggT),(*fNormal)[ii]);
+//              }
+//              if(*ggT>1)//We only need to do this if the ggT is non-trivial
+//              {
+// //                   int64vec *fCopy = fAct->getFacetNormal();
+//                      for(int ii=0;ii<this->numVars;ii++)
+//                              (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT);
+//                      fAct->setFacetNormal(fNormal);
+//              }               
+//              delete fNormal;
+//              delete ggT;
+//              /*And now for the codim2*/
+//              while(codim2Act!=NULL)
+//              {                               
+//                      int64vec *n;
+//                      n=codim2Act->getFacetNormal();
+//                      int *ggT=new int;
+//                      *ggT=1;
+//                      for(int ii=0;ii<this->numVars;ii++)
+//                      {
+//                              *ggT = intgcd((*ggT),(*n)[ii]);
+//                      }
+//                      if(*ggT>1)
+//                      {
+//                              for(int ii=0;ii<this->numVars;ii++)
+//                              {
+//                                      (*n)[ii] = ((*n)[ii])/(*ggT);
+//                              }
+//                              codim2Act->setFacetNormal(n);
+//                      }
+//                      codim2Act = codim2Act->next;
+//                      delete n;
+//                      delete ggT;
+//              }
+//              fAct = fAct->next;
+//      }
+// }
+
+/** \brief Enqueue new facets into the searchlist 
  * The searchlist (SLA for short) is implemented as a FIFO
  * Checks are done as follows:
  * 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
- * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
+ * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the 
  *      facet from SLA and do not add fAct.
  * It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
@@ -3480 +3428 @@
         facet *slHead;
         slHead = f;
-        facet *slAct;   called with f=SearchListRoot
+        facet *slAct;   //called with f=SearchListRoot
         slAct = f;
-        facet *slEnd;   Pointer to end of SLA
+        facet *slEnd;   //Pointer to end of SLA
         slEnd = f;
-        facet *slEndStatic;     //marks the end before a new facet is added
+//      facet *slEndStatic;     //marks the end before a new facet is added             
         facet *fAct;
         fAct = this->facetPtr;
@@ -3494 +3442 @@
         volatile facet *deleteMarker;
         deleteMarker = NULL;
-
+                        
         /** \brief  Flag to mark a facet that might be added
          * The following case may occur:
@@ -3500 +3448 @@
          * This does however not mean that there does not exist a facet behind the current slAct that is actually equal
          * to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
-         * FALSE and the according slAct is deleted.
+         * FALSE and the according slAct is deleted. 
          * If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
          */
 
-        /**A facet was removed, lengthOfSearchlist-- thus we must not rely on
+        /**A facet was removed, lengthOfSearchlist-- thus we must not rely on 
          * if(notParallelCtr==lengthOfSearchList) but rather
          * if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
@@ -3513 +3461 @@
                 slEnd=slEnd->next;
         }
-        /*1st step: compare facetNormals*/
+        /*1st step: compare facetNormals*/                      
         while(fAct!=NULL)
-        {
+        {                                               
                 if(fAct->isFlippable==TRUE)
                 {
                         int64vec *fNormal=NULL;
-                        fNormal=fAct->getFacetNormal();
+                        fNormal=fAct->getFacetNormal();                 
                         slAct = slHead;
-                        /*If slAct==NULL and fAct!=NULL
+                        /*If slAct==NULL and fAct!=NULL 
                         we just copy all remaining facets into SLA*/
                         if(slAct==NULL)
@@ -3529 +3477 @@
                                 while(fCopy!=NULL)
                                 {
-                                        if(fCopy->isFlippable==TRUE)We must assure fCopy is also flippable
+                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
                                         {
                                                 if(slAct==NULL)
                                                 {
-                                                        slAct = new facet(*fCopy/*,TRUE*/);copy constructor
-                                                        slHead = slAct;
+                                                        slAct = new facet(*fCopy/*,TRUE*/);//copy constructor
+                                                        slHead = slAct;                                                         
                                                 }
                                                 else
@@ -3544 +3492 @@
                                         fCopy = fCopy->next;
                                 }
-                                break;Where does this lead to?
+                                break;//Where does this lead to?
                         }
                         /*End of dumping into SLA*/
@@ -3554 +3502 @@
 #ifdef gfan_DEBUG
                                 printf("Checking facet (");fNormal->show(1,1);printf(") against (");slNormal->show(1,1);printf(")\n");
-#endif
-                                if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
-                                        exit(-1);
+#endif                          
+//                              if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
+//                                      exit(-1);
                                 if(areEqual2(fAct,slAct))
-                                {
+                                {                                       
                                         deleteMarker = slAct;
                                         if(slAct==slHead)
-                                        {
-                                                slHead = slAct->next;
+                                        {                                               
+                                                slHead = slAct->next;                                           
                                                 if(slHead!=NULL)
                                                         slHead->prev = NULL;
@@ -3570 +3518 @@
                                                 slEnd=slEnd->prev;
                                                 slEnd->next = NULL;
-                                        }
+                                        }                                                               
                                         else
                                         {
@@ -3580 +3528 @@
                                         if(deleteMarker!=NULL)
                                         {
-                                                delete deleteMarker;
-                                                deleteMarker=NULL;
+//                                              delete deleteMarker;
+//                                              deleteMarker=NULL;
                                         }
 #ifdef gfan_DEBUG
@@ -3587 +3535 @@
 #endif
                                         delete slNormal;
-                                        break;leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
+                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
                                 }
                                 slAct = slAct->next;
@@ -3596 +3544 @@
                                 slAct = slAct->next;*/
                                 if(deleteMarker!=NULL)
-                                {
-                                        delete deleteMarker;
-                                        deleteMarker=NULL;
+                                {                                               
+//                                      delete deleteMarker;
+//                                      deleteMarker=NULL;
                                 }
                                 delete slNormal;
-                                                if slAct was marked as to be deleted, delete it here!
-                        }while(slAct!=NULL)
+                                                //if slAct was marked as to be deleted, delete it here!
+                        }//while(slAct!=NULL)
                         if(removalOccured==FALSE)
                         {
 #ifdef gfan_DEBUG
-                                cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl;
-#endif
-                                Add fAct to SLA
+//                              cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl;
+#endif
+                                //Add fAct to SLA
                                 facet *marker;
                                 marker = slEnd;
@@ -3615 +3563 @@
 
                                 slEnd->next = new facet();
-                                slEnd = slEnd->next;Update slEnd
+                                slEnd = slEnd->next;//Update slEnd
                                 facet *slEndCodim2Root;
                                 facet *slEndCodim2Act;
                                 slEndCodim2Root = slEnd->codim2Ptr;
                                 slEndCodim2Act = slEndCodim2Root;
-
+                                                                
                                 slEnd->setUCN(this->getUCN());
                                 slEnd->isFlippable = TRUE;
                                 slEnd->setFacetNormal(fNormal);
-                                NOTE Add interior point here
-                                This is ugly but needed for flip2
-                                Better: have slEnd->interiorPoint point to fAct->interiorPoint
-                                NOTE Only reference -> c.f. copy constructor
-                                slEnd->setInteriorPoint(fAct->getInteriorPoint());
+                                //NOTE Add interior point here
+                                //This is ugly but needed for flip2
+                                //Better: have slEnd->interiorPoint point to fAct->interiorPoint
+                                //NOTE Only reference -> c.f. copy constructor
+                                //slEnd->setInteriorPoint(fAct->getInteriorPoint());
                                 slEnd->interiorPoint = fAct->interiorPoint;
                                 slEnd->prev = marker;
-                                Copy codim2-facets
-                                int64vec *f2Normal=new int64vec(this->numVars);
+                                //Copy codim2-facets
+                                //int64vec *f2Normal=new int64vec(this->numVars);
                                 while(f2Act!=NULL)
                                 {
@@ -3640 +3588 @@
                                         {
                                                 slEndCodim2Root = new facet(2);
-                                                slEnd->codim2Ptr = slEndCodim2Root;
+                                                slEnd->codim2Ptr = slEndCodim2Root;                     
                                                 slEndCodim2Root->setFacetNormal(f2Normal);
                                                 slEndCodim2Act = slEndCodim2Root;
                                         }
-                                        else
+                                        else                                    
                                         {
                                                 slEndCodim2Act->next = new facet(2);
@@ -3653 +3601 @@
                                         delete f2Normal;
                                 }
-                                gcone::lengthOfSearchList++;
-                        }if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
+                                gcone::lengthOfSearchList++;                                                    
+                        }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
                         fAct = fAct->next;
                         delete fNormal;
-                        delete slNormal;
-                }if(fAct->isFlippable==TRUE)
+//                      delete slNormal;
+                }//if(fAct->isFlippable==TRUE)
                 else
                 {
@@ -3665 +3613 @@
                 if(gcone::maxSize<gcone::lengthOfSearchList)
                         gcone::maxSize= gcone::lengthOfSearchList;
-        }while(fAct!=NULL)
+        }//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
+#endif                                          
         return slHead;
-}addC2N
+}//addC2N
 
 /** Enqueuing using shallow copies*/
@@ -3683 +3631 @@
         facet *slHead;
         slHead = f;
-        facet *slAct;   called with f=SearchListRoot
+        facet *slAct;   //called with f=SearchListRoot
         slAct = f;
-        static facet *slEnd;    Pointer to end of SLA
+        static facet *slEnd;    //Pointer to end of SLA
         if(slEnd==NULL)
                 slEnd = f;
-
+        
         facet *fAct;
-        fAct = this->facetPtr;New facets to compare
+        fAct = this->facetPtr;//New facets to compare
         facet *codim2Act;
         codim2Act = this->facetPtr->codim2Ptr;
@@ -3697 +3645 @@
         {
                 slEnd=slEnd->next;
-        }
+        }       
         while(fAct!=NULL)
         {
@@ -3708 +3656 @@
                                 printf("Zero length SLA\n");
                                 facet *fCopy;
-                                fCopy = fAct;
+                                fCopy = fAct;                           
                                 while(fCopy!=NULL)
                                 {
-                                        if(fCopy->isFlippable==TRUE)We must assure fCopy is also flippable
+                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
                                         {
                                                 if(slAct==NULL)
                                                 {
-                                                        slAct = slAct->shallowCopy(*fCopy);shallow copy constructor
-                                                        slHead = slAct;
+                                                        slAct = slAct->shallowCopy(*fCopy);//shallow copy constructor
+                                                        slHead = slAct;                                                         
                                                 }
                                                 else
@@ -3726 +3674 @@
                                         fCopy = fCopy->next;
                                 }
-                                break;  WTF?
+                                break;  //WTF?
                         }
                         /*Comparison starts here*/
@@ -3734 +3682 @@
 #ifdef gfan_DEBUG
         printf("Checking facet (");fAct->fNormal->show(1,1);printf(") against (");slAct->fNormal->show(1,1);printf(")\n");
-#endif
+#endif  
                                 if(areEqual2(fAct,slAct))
                                 {
@@ -3740 +3688 @@
                                         if(slAct==slHead)
                                         {
-                                                fDeleteMarker=slHead;
-                                                printf("headUCN@enq=%i\n",slHead->getUCN());
-                                                slHead = slAct->next;
-                                                printf("headUCN@enq=%i\n",slHead->getUCN());
+//                                              fDeleteMarker=slHead;
+//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
+                                                slHead = slAct->next;                                           
+//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
                                                 if(slHead!=NULL)
                                                 {
@@ -3749 +3697 @@
                                                 }
                                                 fDeleteMarker->shallowDelete();
-                                                delete fDeleteMarker;//NOTE this messes up fAct in noRevS!
-                                                printf("headUCN@enq=%i\n",slHead->getUCN());
+                                                //delete fDeleteMarker;//NOTE this messes up fAct in noRevS!
+//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
                                         }
                                         else if (slAct==slEnd)
@@ -3758 +3706 @@
                                                 fDeleteMarker->shallowDelete();
                                                 delete(fDeleteMarker);
-                                        }
+                                        }                                                               
                                         else
                                         {
@@ -3771 +3719 @@
 printf("Removing (");fAct->fNormal->show(1,1);printf(") from list\n");
 #endif
-                                        break;leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
+                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
                                 }
-                                slAct = slAct->next;
-                        }while(slAct!=NULL)
+                                slAct = slAct->next;                            
+                        }//while(slAct!=NULL)
                         if(removalOccured==FALSE)
                         {
@@ -3784 +3732 @@
                         }
                         fAct = fAct->next;
-                        if(fDeleteMarker!=NULL)
-                        {
-                                fDeleteMarker->shallowDelete();
-                                delete(fDeleteMarker);
-                                fDeleteMarker=NULL;
-                        }
+//                      if(fDeleteMarker!=NULL)
+//                      {
+//                              fDeleteMarker->shallowDelete();
+//                              delete(fDeleteMarker);
+//                              fDeleteMarker=NULL;
+//                      }
                 }
                 else
                         fAct = fAct->next;
-        }while(fAct!=NULL)
-
+        }//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
-        printf("headUCN@enq=%i\n",slHead->getUCN());
+#endif  
+//      printf("headUCN@enq=%i\n",slHead->getUCN());
         return slHead;
 }
 
-/**
+/** 
 * During flip2 every gc->baseRing gets two ringorder_a
 * To avoid having an awful lot of those in the end we endow
@@ -3824 +3772 @@
         omfree(replacementRing->wvhdl);
         replacementRing->wvhdl =(int **)omAlloc0(4*sizeof(int/*64*/**));
-
+                        
         replacementRing->order[0]=ringorder_a/*64*/;
         replacementRing->block0[0]=1;
         replacementRing->block1[0]=replacementRing->N;
-
+                
         replacementRing->order[1]=ringorder_dp;
         replacementRing->block0[1]=1;
         replacementRing->block1[1]=replacementRing->N;
-
+        
         replacementRing->order[2]=ringorder_C;
 
-        int64vec *ivw = this->getIntPoint(TRUE);returns a reference
-        assert(this->ivIntPt);
-        int length=ivw->length();
+        int64vec *ivw = this->getIntPoint(TRUE);//returns a reference
+//      assert(this->ivIntPt);  
+        int length=ivw->length();       
         int/*64*/ *A=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
         for (int jj=0;jj<length;jj++)
@@ -3843 +3791 @@
                 A[jj]=(*ivw)[jj];
                 if((*ivw)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::replaceDouble_ringorder_a_ByASingleOne!\n");
-        }
-        delete ivw; //Not needed if this->getIntPoint(TRUE)
+        }       
+        //delete ivw; //Not needed if this->getIntPoint(TRUE)
         replacementRing->wvhdl[0]=(int *)A;
         replacementRing->block1[0]=length;
@@ -3853 +3801 @@
         rDelete(this->baseRing);
         this->baseRing=rCopy(replacementRing);
-        this->gcBasis=idCopy(temporaryGroebnerBasis);
+        this->gcBasis=idCopy(temporaryGroebnerBasis);   
         /*And back to where we came from*/
         rChangeCurrRing(srcRing);
         idDelete( (ideal*)&temporaryGroebnerBasis );
-        rDelete(replacementRing);
+        rDelete(replacementRing);       
 }
 
@@ -3877 +3825 @@
         return p;
 }
-
+        
+static int intgcd(const int &a, const int &b)
+{
+        int r, p=a, q=b;
+        if(p < 0)
+                p = -p;
+        if(q < 0)
+                q = -q;
+        while(q != 0)
+        {
+                r = p % q;
+                p = q;
+                q = r;
+        }
+        return p;
+}
+                
 /** \brief Construct a dd_MatrixPtr from a cone's list of facets
  * NO LONGER USED
  */
- inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
- {
-        facet *fAct;
-        fAct = gc.facetPtr;
-
-        dd_MatrixPtr M;
-        dd_rowrange ddrows;
-        dd_colrange ddcols;
-        ddcols=(this->numVars)+1;
-        ddrows=this->numFacets;
-        dd_NumberType numb = dd_Integer;
-        M=dd_CreateMatrix(ddrows,ddcols);
-
-        int jj=0;
-
-        while(fAct!=NULL)
-        {
-                int64vec *comp;
-                comp = fAct->getFacetNormal();
-                for(int ii=0;ii<this->numVars;ii++)
-                {
-                        dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
-                }
-                jj++;
-                delete comp;
-                fAct=fAct->next;
-        }
-        return M;
- }
-
+// inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
+// {
+//      facet *fAct;
+//      fAct = gc.facetPtr;
+//      
+//      dd_MatrixPtr M;
+//      dd_rowrange ddrows;
+//      dd_colrange ddcols;
+//      ddcols=(this->numVars)+1;
+//      ddrows=this->numFacets;
+//      dd_NumberType numb = dd_Integer;
+//      M=dd_CreateMatrix(ddrows,ddcols);                       
+//                      
+//      int jj=0;
+//      
+//      while(fAct!=NULL)
+//      {
+//              int64vec *comp;
+//              comp = fAct->getFacetNormal();
+//              for(int ii=0;ii<this->numVars;ii++)
+//              {
+//                      dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
+//              }
+//              jj++;
+//              delete comp;
+//              fAct=fAct->next;                                
+//      }                       
+//      return M;
+// }
+                
 /** \brief Write information about a cone into a file on disk
  *
@@ -3926 +3890 @@
         stringstream ss;
         ss << UCN;
-        string UCNstr = ss.str();
-
+        string UCNstr = ss.str();               
+                        
         string prefix="/tmp/Singular/cone";
-        string prefix="./";     //crude hack -> run tests in separate directories
+//      string prefix="./";     //crude hack -> run tests in separate directories
         string suffix=".sg";
         string filename;
@@ -3935 +3899 @@
         filename.append(UCNstr);
         filename.append(suffix);
-
-        int thisPID = getpid();
-        ss << thisPID;
-        string strPID = ss.str();
-        string prefix="./";
-
+                
+//      int thisPID = getpid();
+//      ss << thisPID;
+//      string strPID = ss.str();
+//      string prefix="./";
+                                                
         ofstream gcOutputFile(filename.c_str());
         assert(gcOutputFile);
         facet *fAct;
-        fAct = gc.facetPtr;
+        fAct = gc.facetPtr;                     
         facet *f2Act;
         f2Act=fAct->codim2Ptr;
-
+        
         char *ringString = rString(gc.baseRing);
-
+                        
         if (!gcOutputFile)
         {
@@ -3955 +3919 @@
         }
         else
-        {
+        {       
                 gcOutputFile << "UCN" << endl;
                 gcOutputFile << this->UCN << endl;
-                gcOutputFile << "RING" << endl;
+                gcOutputFile << "RING" << endl; 
                 gcOutputFile << ringString << endl;
                 gcOutputFile << "GCBASISLENGTH" << endl;
-                gcOutputFile << IDELEMS(gc.gcBasis) << endl;
-                Write this->gcBasis into file
-                gcOutputFile << "GCBASIS" << endl;
+                gcOutputFile << IDELEMS(gc.gcBasis) << endl;                    
+                //Write this->gcBasis into file
+                gcOutputFile << "GCBASIS" << endl;                              
                 for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
-                {
+                {                                       
                         gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
-                }
-
-                gcOutputFile << "FACETS" << endl;
-
+                }                               
+                                        
+                gcOutputFile << "FACETS" << endl;                                                               
+                
                 while(fAct!=NULL)
-                {
+                {       
                         const int64vec *iv=fAct->getRef2FacetNormal();
-                        iv=fAct->getRef2FacetNormal();//->getFacetNormal();