Changeset 4199b3 in git


Ignore:
Timestamp:
Mar 17, 2011, 5:37:55 PM (13 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
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:
2 edited

Legend:

Unmodified
Added
Removed
  • kernel/gfan.cc

    r9a038bd r4199b3  
    99#include <kernel/mod2.h>
    1010
    11 #ifdef HAVE_GFAN
     11#ifdef HAVE_FANS
    1212#include <kernel/options.h>
    1313#include <kernel/kstd1.h>
    1414#include <kernel/kutil.h>
    15  #include "int64vec.h"
     15// #include "int64vec.h"
    1616#include <kernel/polys.h>
    1717#include <kernel/ideals.h>
    1818#include <kernel/kmatrix.h>
    1919#include <kernel/GMPrat.h>
    20 #include "fast_maps.h"  //Mapping of ideals
    21  #include "maps.h"
    22 #include "ring.h"       apparently not needed
    23  #include "structs.h"
     20
     21#include "ring.h"       //apparently not needed
     22// #include "structs.h"
    2423#include <Singular/lists.h>
    2524#include <kernel/prCopy.h>
    2625#include <kernel/stairc.h>
    27  #include <bitset>
    28 #include <fstream>      read-write cones to files
    29  #include <gmp.h>
     26// #include <bitset>
     27#include <fstream>      //read-write cones to files
     28// #include <gmp.h>
    3029#include <string>
    3130#include <sstream>
    32  #include <time.h>
     31// #include <time.h>
    3332#include <stdlib.h>
    34 #include <gmpxx.h>
    35  #include <vector>
    3633#include <assert.h>
    37 
     34// #include <Singular/bbfan.h>
     35// #include <Singular/bbcone.h>
     36#include <gfanlib/gfanlib.h>
    3837
    3938/*DO NOT REMOVE THIS*/
     
    4241#endif
    4342
    44 Hacks for different working places
     43//Hacks for different working places
    4544#define p800
    46 
    47 #ifdef UNI
    48 #include <ers/urmel/alggeom/monerjan/cddlib/include/setoper.h>
    49 #include <ers/urmel/alggeom/monerjan/cddlib/include/cdd.h>
    50 #endif
    51 
    52 #ifdef HOME
    53 #include <me/momo/studium/diplomarbeit/cddlib/include/setoper.h>
    54 #include <me/momo/studium/diplomarbeit/cddlib/include/cdd.h>
    55 #endif
    56 
    57 #ifdef ITWM
    58 #include <slg/monerjan/cddlib/include/setoper.h>
    59 #include <slg/monerjan/cddlib/include/cdd.h>
    60 #include <slg/monerjan/cddlib/include/cddmp.h>
    61 #endif
    6245
    6346#ifdef p800
     
    6851
    6952#ifndef gfan_DEBUG
    70  #define gfan_DEBUG
     53// #define gfan_DEBUG
    7154#ifndef gfan_DEBUGLEVEL
    7255#define gfan_DEBUGLEVEL 1
     
    7457#endif
    7558
    76 NOTE Defining this will slow things down!
    77 Only good for very coarse profiling
    78  #define gfanp
     59//NOTE Defining this will slow things down!
     60//Only good for very coarse profiling
     61// #define gfanp
    7962#ifdef gfanp
    8063#include <sys/time.h>
     
    8467#ifndef SHALLOW
    8568#define SHALLOW
     69#endif
     70
     71#ifndef USE_ZFAN
     72#define USE_ZFAN
    8673#endif
    8774
     
    9683*
    9784*/
    98 
     85               
    9986/** \brief The default constructor for facets
    10087*/
    101 facet::facet()
    102 {
    103         Pointer to next facet.  */
     88facet::facet()                 
     89{
     90        // Pointer to next facet.  */
    10491        /* Defaults to NULL. This way there is no need to check explicitly */
    10592        this->fNormal=NULL;
    106         this->interiorPoint=NULL;
     93        this->interiorPoint=NULL;               
    10794        this->UCN=0;
    10895        this->codim2Ptr=NULL;
    109         this->codim=1;          default for (codim-1)-facets
     96        this->codim=1;          //default for (codim-1)-facets
    11097        this->numCodim2Facets=0;
    11198        this->numRays=0;
    11299        this->flipGB=NULL;
    113100        this->next=NULL;
    114         this->prev=NULL;
    115         this->flipRing=NULL;    the ring on the other side
     101        this->prev=NULL;               
     102        this->flipRing=NULL;    //the ring on the other side
    116103        this->isFlippable=FALSE;
    117104}
    118 
     105               
    119106/** \brief Constructor for facets of codim == 2
    120107* Note that as of now the code of the constructors is only for facets and codim2-faces. One
     
    124111{
    125112        this->fNormal=NULL;
    126         this->interiorPoint=NULL;
     113        this->interiorPoint=NULL;                       
    127114        this->UCN=0;
    128115        this->codim2Ptr=NULL;
     
    130117        {
    131118                this->codim=n;
    132         }NOTE Handle exception here!
     119        }//NOTE Handle exception here!                 
    133120        this->numCodim2Facets=0;
    134121        this->numRays=0;
     
    139126        this->isFlippable=FALSE;
    140127}
    141 
     128               
    142129/** \brief The copy constructor
    143130* By default only copies the fNormal, f2Normals and UCN
     
    148135        this->UCN=f.UCN;
    149136        this->isFlippable=f.isFlippable;
    150         Needed for flip2
    151         NOTE ONLY REFERENCE
    152         this->interiorPoint=iv64Copy(f.interiorPoint);only referencing is prolly not the best thing to do in a copy constructor
     137        //Needed for flip2
     138        //NOTE ONLY REFERENCE
     139        this->interiorPoint=iv64Copy(f.interiorPoint);//only referencing is prolly not the best thing to do in a copy constructor
    153140        facet* f2Copy;
    154141        f2Copy=f.codim2Ptr;
     
    168155                {
    169156                        f2Act=new facet(2);
    170                         this->codim2Ptr=f2Act;
     157                        this->codim2Ptr=f2Act;                 
    171158                }
    172159                else
     
    180167                int64vec *f2Normal;
    181168                f2Normal = f2Copy->getFacetNormal();
    182                 f2Act->setFacetNormal(f2Copy->getFacetNormal());
     169//              f2Act->setFacetNormal(f2Copy->getFacetNormal());
    183170                f2Act->setFacetNormal(f2Normal);
    184171                delete f2Normal;
    185172                f2Act->setUCN(f2Copy->getUCN());
    186173                f2Copy = f2Copy->next;
    187         }
     174        }       
    188175}
    189176
     
    209196{
    210197#ifdef gfan_DEBUG
    211         printf("shallowdel@UCN %i\n", this->getUCN());
     198//      printf("shallowdel@UCN %i\n", this->getUCN());
    212199#endif
    213200        this->fNormal=NULL;
    214         this->UCN=0;
     201//      this->UCN=0;
    215202        this->interiorPoint=NULL;
    216203        this->codim2Ptr=NULL;
     
    219206        this->flipGB=NULL;
    220207        this->flipRing=NULL;
    221         assert(this->fNormal==NULL);
    222         delete(this);
    223 }
    224 
     208        assert(this->fNormal==NULL);   
     209//      delete(this);
     210}
     211               
    225212/** The default destructor */
    226213facet::~facet()
    227214{
    228215#ifdef gfan_DEBUG
    229         printf("~facet@UCN %i\n",this->getUCN());
     216//      printf("~facet@UCN %i\n",this->getUCN());
    230217#endif
    231218        if(this->fNormal!=NULL)
     
    234221                delete this->interiorPoint;
    235222        /* Cleanup the codim2-structure */
    236         if(this->codim==2)
    237         {
    238                 facet *codim2Ptr;
    239                 codim2Ptr = this->codim2Ptr;
    240                 while(codim2Ptr!=NULL)
    241                 {
    242                         if(codim2Ptr->fNormal!=NULL)
    243                         {
    244                                 delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
    245                                 codim2Ptr = codim2Ptr->next;
    246                         }
    247                 }
    248         }
    249         The rays are stored in the cone!
     223//      if(this->codim==2)
     224//      {
     225//              facet *codim2Ptr;
     226//              codim2Ptr = this->codim2Ptr;
     227//              while(codim2Ptr!=NULL)
     228//              {
     229//                      if(codim2Ptr->fNormal!=NULL)
     230//                      {
     231//                              delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
     232//                              codim2Ptr = codim2Ptr->next;
     233//                      }
     234//              }
     235//      }
     236        //The rays are stored in the cone!
    250237        if(this->flipGB!=NULL)
    251238                idDelete((ideal *)&this->flipGB);
    252         if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
    253                 rDelete(this->flipRing); //See vol II/134
    254         this->flipRing=NULL;
     239//      if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
     240//              rDelete(this->flipRing); //See vol II/134
     241//      this->flipRing=NULL;
    255242        this->prev=NULL;
    256243        this->next=NULL;
    257244}
    258 
     245               
    259246inline const int64vec *facet::getRef2FacetNormal() const
    260247{
    261248        return(this->fNormal);
    262 }
     249}       
    263250
    264251/** Equality check for facets based on unique interior points*/
     
    281268                }
    282269        }
    283         if( fIntP->compare(gIntP)!=0) res=FALSE;
     270//      if( fIntP->compare(gIntP)!=0) res=FALSE;
    284271#ifdef gfanp
    285272        gettimeofday(&end, 0);
    286273        gcone::t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    287 #endif
     274#endif 
    288275        return res;
    289276}
     
    308295        int notParallelCtr=0;
    309296        int ctr=0;
    310         const int64vec* fNormal; No new since iv64Copy and therefore getFacetNormal return a new
     297        const int64vec* fNormal; //No new since iv64Copy and therefore getFacetNormal return a new
    311298        const int64vec* sNormal;
    312         fNormal = f->getRef2FacetNormal();->getFacetNormal();
    313         sNormal = s->getRef2FacetNormal();->getFacetNormal();
    314         Do not need parallelity. Too time consuming
    315         if(!isParallel(*fNormal,*sNormal))
    316         if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
    317                 notParallelCtr++;
    318         else//parallelity, so we check the codim2-facets
     299        fNormal = f->getRef2FacetNormal();//->getFacetNormal();
     300        sNormal = s->getRef2FacetNormal();//->getFacetNormal();
     301#include "intvec.h"
     302        //Do not need parallelity. Too time consuming
     303//      if(!isParallel(*fNormal,*sNormal))
     304//      if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
     305 //             notParallelCtr++;
     306//      else//parallelity, so we check the codim2-facets
    319307        int64vec *fNRef=const_cast<int64vec*>(fNormal);
    320308        int64vec *sNRef=const_cast<int64vec*>(sNormal);
    321         if(isParallel(*fNormal,*sNormal))
     309//      if(isParallel(*fNormal,*sNormal))       
    322310        if(isParallel(*fNRef,*sNRef))
    323         if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
     311//      if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
    324312        {
    325313                facet* f2Act;
    326314                facet* s2Act;
    327                 f2Act = f->codim2Ptr;
     315                f2Act = f->codim2Ptr;           
    328316                ctr=0;
    329317                while(f2Act!=NULL)
    330318                {
    331319                        const int64vec* f2Normal;
    332                         f2Normal = f2Act->getRef2FacetNormal();->getFacetNormal();
    333                         int64vec *f2Ref=const_cast<int64vec*>(f2Normal);
     320                        f2Normal = f2Act->getRef2FacetNormal();//->getFacetNormal();
     321//                      int64vec *f2Ref=const_cast<int64vec*>(f2Normal);
    334322                        s2Act = s->codim2Ptr;
    335323                        while(s2Act!=NULL)
    336324                        {
    337325                                const int64vec* s2Normal;
    338                                 s2Normal = s2Act->getRef2FacetNormal();->getFacetNormal();
    339                                 bool foo=areEqual(f2Normal,s2Normal);
    340                                 int64vec *s2Ref=const_cast<int64vec*>(s2Normal);
     326                                s2Normal = s2Act->getRef2FacetNormal();//->getFacetNormal();
     327//                              bool foo=areEqual(f2Normal,s2Normal);
     328//                              int64vec *s2Ref=const_cast<int64vec*>(s2Normal);
    341329                                int foo=f2Normal->compare(s2Normal);
    342330                                if(foo==0)
    343331                                        ctr++;
    344332                                s2Act = s2Act->next;
    345                                 delete s2Normal;
     333//                              delete s2Normal;
    346334                        }
    347                         delete f2Normal;
     335//                      delete f2Normal;
    348336                        f2Act = f2Act->next;
    349                 }
    350         }
    351         delete fNormal;
    352         delete sNormal;
     337                }               
     338        }
     339//      delete fNormal;
     340//      delete sNormal;
    353341        if(ctr==f->numCodim2Facets)
    354342                res=TRUE;
     
    365353#endif
    366354        return res;
    367         int64vec *foo=ivNeg(sNormal);
    368         if(fNormal->compare(foo)!=0) //facet normals
    369         {
    370                 delete foo;
    371                 res=FALSE;
    372         }
    373         else
    374         {
    375                 facet* f2Act;
    376                 facet* s2Act;
    377                 f2Act = f->codim2Ptr;
    378                 ctr=0;
    379                 while(f2Act!=NULL)
    380                 {
    381                         int64vec* f2Normal;
    382                         f2Normal = f2Act->getFacetNormal();
    383                         s2Act = s->codim2Ptr;
    384                         while(s2Act!=NULL)
    385                         {
    386                                 int64vec* s2Normal;
    387                                 s2Normal = s2Act->getFacetNormal();
    388  //                             bool foo=areEqual(f2Normal,s2Normal);
    389                                 int foo=f2Normal->compare(s2Normal);
    390                                 if(foo==0)
    391                                         ctr++;
    392                                 s2Act = s2Act->next;
    393                                 delete s2Normal;
    394                         }
    395                         delete f2Normal;
    396                         f2Act = f2Act->next;
    397                 }
    398         }
    399         delete fNormal;
    400         delete sNormal;
    401         if(ctr==f->numCodim2Facets)
    402                 res=TRUE;
    403         return res;
    404 }
    405 
     355//      int64vec *foo=ivNeg(sNormal);
     356//      if(fNormal->compare(foo)!=0) //facet normals
     357//      {
     358//              delete foo;
     359//              res=FALSE;
     360//      }
     361//      else
     362//      {
     363//              facet* f2Act;
     364//              facet* s2Act;
     365//              f2Act = f->codim2Ptr;           
     366//              ctr=0;
     367//              while(f2Act!=NULL)
     368//              {
     369//                      int64vec* f2Normal;
     370//                      f2Normal = f2Act->getFacetNormal();
     371//                      s2Act = s->codim2Ptr;
     372//                      while(s2Act!=NULL)
     373//                      {
     374//                              int64vec* s2Normal;
     375//                              s2Normal = s2Act->getFacetNormal();
     376// //                           bool foo=areEqual(f2Normal,s2Normal);
     377//                              int foo=f2Normal->compare(s2Normal);
     378//                              if(foo==0)
     379//                                      ctr++;
     380//                              s2Act = s2Act->next;
     381//                              delete s2Normal;
     382//                      }
     383//                      delete f2Normal;
     384//                      f2Act = f2Act->next;
     385//              }
     386//      }               
     387//      delete fNormal;
     388//      delete sNormal;
     389//      if(ctr==f->numCodim2Facets)
     390//              res=TRUE;
     391//      return res;
     392}       
     393               
    406394/** Stores the facet normal \param int64vec*/
    407395inline void facet::setFacetNormal(int64vec *iv)
     
    409397        if(this->fNormal!=NULL)
    410398                delete this->fNormal;
    411         this->fNormal = iv64Copy(iv);
    412 }
    413 
    414 /** Hopefully returns the facet normal
     399        this->fNormal = iv64Copy(iv);                   
     400}
     401               
     402/** Hopefully returns the facet normal 
    415403* Mind: iv64Copy returns a new int64vec, so use this in the following way:
    416404* int64vec *iv;
     
    420408*/
    421409inline int64vec *facet::getFacetNormal() const
    422 {
     410{                               
    423411        return iv64Copy(this->fNormal);
    424         return this->fNormal;
     412//      return this->fNormal;
    425413}
    426414
     
    430418        fNormal->show();
    431419}
    432 
     420               
    433421/** Store the flipped GB*/
    434422inline void facet::setFlipGB(ideal I)
     
    436424        this->flipGB=idCopy(I);
    437425}
    438 
     426               
    439427/** Returns a pointer to the flipped GB
    440428Seems not be used
     
    445433        return this->flipGB;
    446434}
    447 
     435               
    448436/** Print the flipped GB*/
    449437inline void facet::printFlipGB()
     
    453441#endif
    454442}
    455 
     443               
    456444/** Set the UCN */
    457445inline void facet::setUCN(int n)
     
    459447        this->UCN=n;
    460448}
    461 
    462 /** \brief Get the UCN
     449               
     450/** \brief Get the UCN 
    463451* Returns the UCN iff this != NULL, else -1
    464452*/
     
    474462#ifdef NDEBUG
    475463        if(this!=NULL)
    476 #endif
     464#endif                 
    477465                return this->UCN;
    478466        else
    479467                return -1;
    480468}
    481 
     469               
    482470/** Store an interior point of the facet */
    483471inline void facet::setInteriorPoint(int64vec *iv)
     
    487475        this->interiorPoint = iv64Copy(iv);
    488476}
    489 
     477               
    490478/** Returns a pointer to this->interiorPoint
    491479* MIND: iv64Copy returns a new int64vec
     
    501489        return (this->interiorPoint);
    502490}
    503 
     491               
    504492/** \brief Debugging function
    505493* prints the facet normal an all (codim-2)-facets that belong to it
    506494*/
    507495volatile void facet::fDebugPrint()
    508 {
    509         facet *codim2Act;
     496{                       
     497        facet *codim2Act;                       
    510498        codim2Act = this->codim2Ptr;
    511499        int64vec *fNormal;
    512         fNormal = this->getFacetNormal();
     500        fNormal = this->getFacetNormal();       
    513501        printf("=======================\n");
    514502        printf("Facet normal = (");fNormal->show(1,1);printf(")\n");
     
    524512        }
    525513        printf("=======================\n");
    526         delete fNormal;
    527 }
    528 
    529 friend class gcone;     //Bad style
     514        delete fNormal; 
     515}               
     516               
     517//friend class gcone;   //Bad style
    530518
    531519
     
    538526
    539527
    540 /** \brief Default constructor.
     528/** \brief Default constructor. 
    541529 *
    542530 * Initialises this->next=NULL and this->facetPtr=NULL
     
    546534        this->next=NULL;
    547535        this->prev=NULL;
    548         this->facetPtr=NULL;    maybe this->facetPtr = new facet();
     536        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
    549537        this->baseRing=currRing;
    550538        this->counter++;
    551         this->UCN=this->counter;
     539        this->UCN=this->counter;                       
    552540        this->numFacets=0;
    553541        this->ivIntPt=NULL;
    554542        this->gcRays=NULL;
    555543}
    556 
     544               
    557545/** \brief Constructor with ring and ideal
    558546 *
    559  * This constructor takes the root ring and the root ideal as parameters and stores
     547 * This constructor takes the root ring and the root ideal as parameters and stores 
    560548 * them in the private members gcone::rootRing and gcone::inputIdeal
    561549 * This constructor is only called once in the computation of the Gröbner fan,
     
    569557        this->next=NULL;
    570558        this->prev=NULL;
    571         this->facetPtr=NULL;
     559        this->facetPtr=NULL;                   
    572560        this->inputIdeal=I;
    573561        this->baseRing=currRing;
     
    579567        this->gcRays=NULL;
    580568}
    581 
    582 /** \brief Copy constructor
     569               
     570/** \brief Copy constructor 
    583571 *
    584572 * Copies a cone, sets this->gcBasis to the flipped GB
    585573 * Call this only after a successful call to gcone::flip which sets facet::flipGB
    586 */
     574*/             
    587575gcone::gcone(const gcone& gc, const facet &f)
    588576{
    589577        this->next=NULL;
    590         this->prev=(gcone *)&gc; //comment in to get a tree
     578//      this->prev=(gcone *)&gc; //comment in to get a tree
    591579        this->prev=NULL;
    592         this->numVars=gc.numVars;
     580        this->numVars=gc.numVars;                                               
    593581        this->counter++;
    594582        this->UCN=this->counter;
     
    596584        this->facetPtr=NULL;
    597585        this->gcBasis=idrCopyR(f.flipGB, f.flipRing);
    598         this->inputIdeal=idCopy(this->gcBasis);
     586//      this->inputIdeal=idCopy(this->gcBasis);
    599587        this->baseRing=rCopy(f.flipRing);
    600588        this->numFacets=0;
     
    602590        this->gcRays=NULL;
    603591}
    604 
    605 /** \brief Default destructor
     592               
     593/** \brief Default destructor 
    606594*/
    607595gcone::~gcone()
     
    619607                idDelete((ideal *)&this->gcBasis);
    620608#endif
    621         idDelete((ideal *)&this->gcBasis);
    622         if(this->inputIdeal!=NULL)
    623                 idDelete((ideal *)&this->inputIdeal);
    624         if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
    625                 rDelete(this->rootRing);
     609//      idDelete((ideal *)&this->gcBasis);
     610//      if(this->inputIdeal!=NULL)
     611//              idDelete((ideal *)&this->inputIdeal);
     612//      if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
     613//              rDelete(this->rootRing);
    626614        if(this->UCN!=1 && this->baseRing!=NULL)
    627615                rDelete(this->baseRing);
     
    638626        }
    639627        this->counter--;
    640         should be deleted in noRevS
    641         dd_FreeMatrix(this->ddFacets);
    642         dd_FreeMatrix(this->ddFacets);
     628        //should be deleted in noRevS
     629//      dd_FreeMatrix(this->ddFacets);
     630        //dd_FreeMatrix(this->ddFacets);
    643631        for(int ii=0;ii<this->numRays;ii++)
    644632                delete(gcRays[ii]);
    645633        omFree(gcRays);
    646 }
     634}                       
    647635
    648636/** Returns the number of cones existing at the time*/
     
    651639        return this->counter;
    652640}
    653 
     641       
    654642/** \brief Set the interior point of a cone */
    655643inline void gcone::setIntPoint(int64vec *iv)
     
    659647        this->ivIntPt=iv64Copy(iv);
    660648}
    661 
     649               
    662650/** \brief Returns either a physical copy the interior point of a cone or just a reference to it.*/
    663651inline int64vec *gcone::getIntPoint(bool shallow)
     
    668656                return iv64Copy(this->ivIntPt);
    669657}
    670 
     658               
    671659/** \brief Print the interior point */
    672660inline void gcone::showIntPoint()
     
    674662        ivIntPt->show();
    675663}
    676 
     664               
    677665/** \brief Print facets
    678666 * This is mainly for debugging purposes. Usually called from within gdb
     
    708696        printf("\n");
    709697}
    710 
     698               
    711699/** For debugging purposes only */
    712700static volatile void showSLA(facet &f)
     
    718706                facet *codim2Act;
    719707                codim2Act = fAct->codim2Ptr;
    720 
     708               
    721709                printf("\n");
    722710                while(fAct!=NULL)
    723711                {
    724                         int64vec *fNormal;
     712                        int64vec *fNormal;             
    725713                        fNormal=fAct->getFacetNormal();
    726714                        printf("(");fNormal->show(1,0);
     
    745733        }
    746734}
    747 
     735               
    748736static void idDebugPrint(const ideal &I)
    749737{
     
    761749static void invPrint(const ideal &I)
    762750{
    763         int numElts=IDELEMS(I);
    764         cout << "inv = ";
    765         for(int ii=0;ii<numElts;ii++);
    766         {
    767                 pWrite0(pHead(I->m[ii]));
    768                 cout << ",";
    769         }
    770         cout << endl;
     751//      int numElts=IDELEMS(I);
     752//      cout << "inv = ";
     753//      for(int ii=0;ii<numElts;ii++);
     754//      {
     755//              pWrite0(pHead(I->m[ii]));
     756//              cout << ",";
     757//      }
     758//      cout << endl;
    771759}
    772760
     
    784772        return res;
    785773}
    786 
     774               
    787775/** \brief Set gcone::numFacets */
    788776inline void gcone::setNumFacets()
    789777{
    790778}
    791 
     779               
    792780/** \brief Get gcone::numFacets */
    793781inline int gcone::getNumFacets()
     
    795783        return this->numFacets;
    796784}
    797 
     785               
    798786inline int gcone::getUCN()
    799787{
    800         if( this!=NULL) && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
     788        if( this!=NULL)// && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
    801789                return this->UCN;
    802790        else
     
    829817 * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
    830818 * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
    831  * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
     819 * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects 
    832820 * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
    833821 * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
     
    844832#endif
    845833        poly aktpoly;
    846         int rows;                        will contain the dimensions of the ineq matrix - deprecated by
     834        int rows;                       // will contain the dimensions of the ineq matrix - deprecated by
    847835        dd_rowrange ddrows;
    848836        dd_colrange ddcols;
    849         dd_rowset ddredrows;            # of redundant rows in ddineq
    850         dd_rowset ddlinset;              the opposite
    851         dd_rowindex ddnewpos=NULL;                all to make dd_Canonicalize happy
    852         dd_NumberType ddnumb=dd_Integer; Number type
    853         dd_ErrorType dderr=dd_NoError;
    854         Compute the # inequalities i.e. rows of the matrix
    855         rows=0; Initialization
     837        dd_rowset ddredrows;            // # of redundant rows in ddineq
     838        dd_rowset ddlinset;             // the opposite
     839        dd_rowindex ddnewpos=NULL;                // all to make dd_Canonicalize happy
     840        dd_NumberType ddnumb=dd_Integer; //Number type
     841        dd_ErrorType dderr=dd_NoError;         
     842        //Compute the # inequalities i.e. rows of the matrix
     843        rows=0; //Initialization
    856844        for (int ii=0;ii<IDELEMS(I);ii++)
    857845        {
    858                 aktpoly=(poly)I->m[ii];
    859                 rows=rows+pLength(aktpoly)-1;
     846//              aktpoly=(poly)I->m[ii];
     847//              rows=rows+pLength(aktpoly)-1;
    860848                rows=rows+pLength((poly)I->m[ii])-1;
    861849        }
    862850
    863         dd_rowrange aktmatrixrow=0;      needed to store the diffs of the expvects in the rows of ddineq
     851        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
    864852        ddrows=rows;
    865853        ddcols=this->numVars;
    866         dd_MatrixPtr ddineq;            Matrix to store the inequalities
    867         ddineq=dd_CreateMatrix(ddrows,ddcols+1); The first col has to be 0 since cddlib checks for additive consts there
    868 
    869         We loop through each g\in GB and compute the resulting inequalities
     854        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
     855        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
     856               
     857        // We loop through each g\in GB and compute the resulting inequalities
    870858        for (int i=0; i<IDELEMS(I); i++)
    871859        {
    872                 aktpoly=(poly)I->m[i];          get aktpoly as i-th component of I
    873                 simpler version of storing expvect diffs
     860                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
     861                //simpler version of storing expvect diffs
    874862                int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
    875863                pGetExpV(aktpoly,leadexpv);
     
    879867                        pNextTerm/*aktpoly*/=pNext(pNextTerm);
    880868                        int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
    881                         pGetExpV(pNextTerm,tailexpv);
     869                        pGetExpV(pNextTerm,tailexpv);                   
    882870                        for(int kk=1;kk<=this->numVars;kk++)
    883                         {
     871                        {                               
    884872                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
    885873                        }
    886874                        aktmatrixrow += 1;
    887875                        omFree(tailexpv);
    888                 }
    889                 omFree(leadexpv);
    890         } for
     876                }               
     877                omFree(leadexpv);       
     878        } //for
    891879#if true
    892880        /*Let's make the preprocessing here. This could already be done in the above for-loop,
     
    895883        * Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
    896884        */
    897         ideal initialForm=idInit(IDELEMS(I),1);
    898         int64vec *gamma=new int64vec(this->numVars);
     885//      ideal initialForm=idInit(IDELEMS(I),1);
     886//      int64vec *gamma=new int64vec(this->numVars);
    899887        int falseGammaCounter=0;
    900888        int *redRowsArray=NULL;
    901889        int num_alloc=0;
    902         int num_elts=0;
     890        int num_elts=0; 
    903891        for(int ii=0;ii<ddineq->rowsize;ii++)
    904892        {
    905893                ideal initialForm=idInit(IDELEMS(I),I->rank);
    906                 read row ii into gamma
    907                 int64 tmp;
     894                //read row ii into gamma
     895//              int64 tmp;
    908896                int64vec *gamma=new int64vec(this->numVars);
    909897                for(int jj=1;jj<=this->numVars;jj++)
     
    915903                computeInv((ideal&)I,initialForm,*gamma);
    916904                delete gamma;
    917                 Create leading ideal
     905                //Create leading ideal
    918906                ideal L=idInit(IDELEMS(initialForm),1);
    919907                for(int jj=0;jj<IDELEMS(initialForm);jj++)
     
    922910                        L->m[jj]=pCopy(/*pHead(initialForm->m[jj]))*/p);
    923911                        pDelete(&p);
    924                 }
    925 
     912                }               
     913               
    926914                LObject *P = new sLObject();//TODO What's the difference between sLObject and LObject?
    927915                memset(P,0,sizeof(LObject));
     
    930918                {
    931919                        bool isMaybeFacet=FALSE;
    932                         P->p1=initialForm->m[jj];       build spolys of initialForm in_v
     920                        P->p1=initialForm->m[jj];       //build spolys of initialForm in_v
    933921
    934922                        for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
    935923                        {
    936924                                P->p2=initialForm->m[kk];
    937                                 ksCreateSpoly(P);
    938                                 if(P->p!=NULL)  spoly non zero=?
    939                                 {
     925                                ksCreateSpoly(P);                                                       
     926                                if(P->p!=NULL)  //spoly non zero=?
     927                                {       
    940928                                        poly p;//NOTE Don't use pInit here. Evil memleak will follow
    941929                                        poly q;
    942930                                        poly pDel,qDel;
    943931                                        p=pCopy(P->p);
    944                                         q=pHead(p);     Monomial q
     932                                        q=pHead(p);     //Monomial q
    945933                                        pDelete(&q);
    946934                                        pDel=p; qDel=q;
    947935                                        isMaybeFacet=FALSE;
    948                                         TODO: Suffices to check LTs here
     936                                        //TODO: Suffices to check LTs here
    949937                                        while(p!=NULL)
    950                                         {
     938                                        {                                               
    951939                                                q=pHead(p);
    952940                                                for(int ll=0;ll<IDELEMS(L);ll++)
    953941                                                {
    954942                                                        if(pLmEqual(L->m[ll],q) || pDivisibleBy(L->m[ll],q))
    955                                                         {
     943                                                        {                                                       
    956944                                                                isMaybeFacet=TRUE;
    957                                                                 break;for
     945                                                                break;//for
    958946                                                        }
    959947                                                }
     
    961949                                                if(isMaybeFacet==TRUE)
    962950                                                {
    963                                                         break;while(p!=NULL)
     951                                                        break;//while(p!=NULL)
    964952                                                }
    965953                                                p=pNext(p);
    966                                         }while
    967                                         pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
     954                                        }//while
     955//                                      pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
    968956                                        if(q!=NULL) pDelete(&q);
    969957                                        pDelete(&pDel);
     
    971959                                        if(isMaybeFacet==FALSE)
    972960                                        {
    973                                                 dd_set_si(ddineq->matrix[ii][0],1);
    974                                                 if(num_alloc==0)
    975                                                         num_alloc += 1;
    976                                                 else
    977                                                         num_alloc += 1;
     961                                                dd_set_si(ddineq->matrix[ii][0],1);                                             
     962//                                              if(num_alloc==0)
     963//                                                      num_alloc += 1;
     964//                                              else                                           
     965//                                                      num_alloc += 1;
    978966                                                if(num_alloc==num_elts) num_alloc==0 ? num_alloc=1 : num_alloc*=2;
    979 
     967                                               
    980968                                                void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
    981969                                                if(!tmp)
     
    987975                                                redRowsArray[num_elts]=ii;
    988976                                                num_elts++;
    989                                                 break;//for(int kk, since we have found one that is not in L
    990                                                 goto _start;    mea culpa, mea culpa, mea maxima culpa
     977                                                //break;//for(int kk, since we have found one that is not in L 
     978                                                goto _start;    //mea culpa, mea culpa, mea maxima culpa
    991979                                        }
    992                                 }if(P->p!=NULL)
     980                                }//if(P->p!=NULL)
    993981                                pDelete(&(P->p));
    994                         }for k
    995                 }for jj
     982                        }//for k
     983                }//for jj
    996984                _start:;
    997985                idDelete(&L);
    998986                delete P;
    999987                idDelete(&initialForm);
    1000         }for(ii<ddineq-rowsize
    1001         delete gamma;
    1002         int offset=0;needed for correction of redRowsArray[ii]
     988        }//for(ii<ddineq-rowsize
     989//      delete gamma;
     990        int offset=0;//needed for correction of redRowsArray[ii]
    1003991#ifdef gfan_DEBUG
    1004992        printf("Removed %i of %i in preprocessing step\n",num_elts,ddineq->rowsize);
     
    1006994        for( int ii=0;ii<num_elts;ii++ )
    1007995        {
    1008                 dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);cddlib sucks at enumeration
     996                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
    1009997                offset++;
    1010         }
    1011         free(redRowsArray);NOTE May crash on some machines.
     998        }       
     999        free(redRowsArray);//NOTE May crash on some machines.
    10121000        /*And now for the strictly positive rows
    10131001        * Doesn't gain significant speedup
     
    10221010                        (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
    10231011                bool isStrictlyPos=FALSE;
    1024                 int posCtr=0;
     1012                int posCtr=0;           
    10251013                for(int jj=0;jj<this->numVars;jj++)
    10261014                {
     
    10301018                        {
    10311019                                if ((*ivPos)[jj]>=0)
    1032                                 {
    1033                                         posCtr++;
     1020                                {                               
     1021                                        posCtr++;                               
    10341022                                }
    1035                         }
     1023                        }                       
    10361024                        delete ivCanonical;
    10371025                }
     
    10521040                        posRowsArray = (int*)tmp;
    10531041                        posRowsArray[num_elts]=ii;
    1054                         num_elts++;
     1042                        num_elts++;     
    10551043                }
    10561044                delete ivPos;
     
    10651053#endif
    10661054
    1067         dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
    1068         ddrows = ddineq->rowsize;       Size of the matrix with redundancies removed
     1055        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);       
     1056        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
    10691057        ddcols = ddineq->colsize;
    1070 
     1058       
    10711059        this->ddFacets = dd_CopyMatrix(ddineq);
    1072 
    1073         /*Write the normals into class facet*/
    1074         facet *fAct;    pointer to active facet
     1060                       
     1061        /*Write the normals into class facet*/ 
     1062        facet *fAct;    //pointer to active facet       
    10751063        int numNonFlip=0;
    10761064        for (int kk = 0; kk<ddrows; kk++)
    10771065        {
    1078                 int64 ggT=1;NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
    1079                 int64vec *load = new int64vec(this->numVars);int64vec to store a single facet normal that will then be stored via setFacetNormal
     1066                int64 ggT=1;//NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
     1067                int64vec *load = new int64vec(this->numVars);//int64vec to store a single facet normal that will then be stored via setFacetNormal
    10801068                for (int jj = 1; jj <ddcols; jj++)
    10811069                {
    10821070                        int64 val;
    10831071                        val = (int64)mpq_get_d(ddineq->matrix[kk][jj]);
    1084                         (*load)[jj-1] = val;    store typecasted entry at pos jj-1 of load
     1072                        (*load)[jj-1] = val;    //store typecasted entry at pos jj-1 of load
    10851073                        ggT = int64gcd(ggT,/*(int64&)foo*/val);
    1086                 }for (int jj = 1; jj <ddcols; jj++)
     1074                }//for (int jj = 1; jj <ddcols; jj++)
    10871075                if(ggT>1)
    10881076                {
    10891077                        for(int ll=0;ll<this->numVars;ll++)
    1090                                 (*load)[ll] /= ggT;make primitive vector
     1078                                (*load)[ll] /= ggT;//make primitive vector
    10911079                }
    10921080                /*Quick'n'dirty hack for flippability. Executed only if gcone::hasHomInput==FALSE
    10931081                * Otherwise every facet intersects the positive orthant
    1094                 */
     1082                */     
    10951083                if(gcone::hasHomInput==FALSE)
    10961084                {
    1097                         TODO: No dP needed
     1085                        //TODO: No dP needed
    10981086                        bool isFlip=FALSE;
    10991087                        for(int jj = 0; jj<load->length(); jj++)
    11001088                        {
    1101                                 int64vec *ivCanonical = new int64vec(load->length());
    1102                                 (*ivCanonical)[jj]=1;
    1103                                 if (dotProduct(*load,*ivCanonical)<0)
    1104                                 {
    1105                                         isFlip=TRUE;
    1106                                         break;  //URGHS
    1107                                 }
    1108                                 delete ivCanonical;
     1089//                              int64vec *ivCanonical = new int64vec(load->length());
     1090//                              (*ivCanonical)[jj]=1;
     1091//                              if (dotProduct(*load,*ivCanonical)<0)   
     1092//                              {
     1093//                                      isFlip=TRUE;
     1094//                                      break;  //URGHS
     1095//                              }
     1096//                              delete ivCanonical;
    11091097                                if((*load)[jj]<0)
    11101098                                {
    11111099                                        isFlip=TRUE;
    11121100                                        break;
    1113                                 }
     1101                                }                               
    11141102                        }/*End of check for flippability*/
    1115                         if(iv64isStrictlyPositive(load))
    1116                                 isFlip=TRUE;
     1103//                      if(iv64isStrictlyPositive(load))
     1104//                              isFlip=TRUE;
    11171105                        if(isFlip==FALSE)
    11181106                        {
     
    11231111                                        facet *fRoot = new facet();
    11241112                                        this->facetPtr = fRoot;
    1125                                         fAct = fRoot;
     1113                                        fAct = fRoot;                                                   
    11261114                                }
    11271115                                else
     
    11341122                                fAct->setUCN(this->getUCN());
    11351123#ifdef gfan_DEBUG
    1136                                 printf("Marking facet (");load->show(1,0);printf(") as non flippable\n");
     1124                                printf("Marking facet (");load->show(1,0);printf(") as non flippable\n");               
    11371125#endif
    11381126                        }
     
    11531141                                fAct->isFlippable=TRUE;
    11541142                                fAct->setFacetNormal(load);
    1155                                 fAct->setUCN(this->getUCN());
     1143                                fAct->setUCN(this->getUCN());                                   
    11561144                        }
    1157                 }hasHomInput==FALSE
    1158                 else    Every facet is flippable
    1159                 {       /*Now load should be full and we can call setFacetNormal*/
     1145                }//hasHomInput==FALSE
     1146                else    //Every facet is flippable
     1147                {       /*Now load should be full and we can call setFacetNormal*/                                     
    11601148                        this->numFacets++;
    11611149                        if(this->numFacets==1)
     
    11721160                        fAct->isFlippable=TRUE;
    11731161                        fAct->setFacetNormal(load);
    1174                         fAct->setUCN(this->getUCN());
    1175                 }if (isFlippable==FALSE)
    1176                 delete load;
    1177         }for (int kk = 0; kk<ddrows; kk++)
    1178 
    1179         In cases like I=<x-1,y-1> there are only non-flippable facets...
     1162                        fAct->setUCN(this->getUCN());                                   
     1163                }//if (isFlippable==FALSE)
     1164                delete load;                           
     1165        }//for (int kk = 0; kk<ddrows; kk++)
     1166                       
     1167        //In cases like I=<x-1,y-1> there are only non-flippable facets...
    11801168        if(numNonFlip==this->numFacets)
    1181         {
     1169        {                                       
    11821170                WerrorS ("Only non-flippable facets. Terminating...\n");
    1183                 exit(-1);//Bit harsh maybe...
    1184         }
    1185 
     1171//              exit(-1);//Bit harsh maybe...
     1172        }
     1173                       
    11861174        /*
    11871175        Now we should have a linked list containing the facet normals of those facets that are
     
    11891177        -flipable
    11901178        Adressing is done via *facetPtr
    1191         */
     1179        */                     
    11921180        if (compIntPoint==TRUE)
    11931181        {
     
    11981186                {
    11991187                        dd_set_si(posRestr->matrix[ii][jj],1);
    1200                         jj++;
     1188                        jj++;                                                   
    12011189                }
    12021190                dd_MatrixAppendTo(&ddineq,posRestr);
    1203                 interiorPoint(ddineq, *iv);     NOTE ddineq contains non-flippable facets
    1204                 this->setIntPoint(iv);  stores the interior point in gcone::ivIntPt
     1191                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
     1192                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
    12051193                delete iv;
    12061194                dd_FreeMatrix(posRestr);
    1207         }
    1208         //Clean up but don't delete the return value!
    1209         dd_FreeMatrix(ddineq);
    1210         set_free(ddredrows);check
    1211         set_free(ddlinset);check
    1212         free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
     1195        }       
     1196        //Clean up but don't delete the return value!   
     1197        //dd_FreeMatrix(ddineq);       
     1198        set_free(ddredrows);//check
     1199        set_free(ddlinset);//check
     1200        //free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
    12131201#ifdef gfanp
    12141202        gettimeofday(&end, 0);
     
    12161204#endif
    12171205
    1218 }gcone::getConeNormals(ideal I)
    1219 
     1206}//gcone::getConeNormals(ideal I)
     1207               
    12201208/** \brief Compute the (codim-2)-facets of a given cone
    12211209 * This method is used during noRevS
     
    12291217        gettimeofday(&start, 0);
    12301218#endif
    1231         this->facetPtr->codim2Ptr = new facet(2);       //instantiate a (codim-2)-facet
     1219        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
    12321220        facet *fAct;
    1233         fAct = this->facetPtr;
     1221        fAct = this->facetPtr;         
    12341222        facet *codim2Act;
    1235         codim2Act = this->facetPtr->codim2Ptr;
    1236         dd_MatrixPtr ddineq;,P,ddakt;
     1223        //codim2Act = this->facetPtr->codim2Ptr;
     1224        dd_MatrixPtr ddineq;//,P,ddakt;
    12371225        dd_ErrorType err;
    1238         ddineq = facets2Matrix(gc);     //get a matrix representation of the cone
     1226        //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
    12391227        ddineq = dd_CopyMatrix(gc.ddFacets);
    12401228        /*Now set appropriate linearity*/
    1241         for (int ii=0; ii<this->numFacets; ii++)
    1242         {
     1229        for (int ii=0; ii<this->numFacets; ii++)                       
     1230        {       
    12431231                dd_rowset impl_linset, redset;
    12441232                dd_rowindex newpos;
    12451233                dd_MatrixPtr ddakt;
    12461234                ddakt = dd_CopyMatrix(ddineq);
    1247                 ddakt->representation=dd_Inequality;    //Not using this makes it faster. But why does the quick check below still work?
    1248                 ddakt->representation=dd_Generator;
     1235//              ddakt->representation=dd_Inequality;    //Not using this makes it faster. But why does the quick check below still work?
     1236//              ddakt->representation=dd_Generator;
    12491237                set_addelem(ddakt->linset,ii+1);/*Now set appropriate linearity*/
    12501238#ifdef gfanp
    12511239                timeval t_ddMC_start, t_ddMC_end;
    12521240                gettimeofday(&t_ddMC_start,0);
    1253 #endif
    1254                 dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
     1241#endif                         
     1242                //dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
    12551243                dd_PolyhedraPtr ddpolyh;
    12561244                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    1257                 ddpolyh=dd_DDMatrix2Poly2(ddakt, dd_MaxCutoff, &err);
     1245//              ddpolyh=dd_DDMatrix2Poly2(ddakt, dd_MaxCutoff, &err);
    12581246                dd_MatrixPtr P;
    1259                 P=dd_CopyGenerators(ddpolyh);
     1247                P=dd_CopyGenerators(ddpolyh);           
    12601248                dd_FreePolyhedra(ddpolyh);
    1261                 TODO Call for one cone , normalize - check equalities - plus lineality -done
     1249                //TODO Call for one cone , normalize - check equalities - plus lineality -done
    12621250#ifdef gfanp
    12631251                gettimeofday(&t_ddMC_end,0);
    12641252                t_ddMC += (t_ddMC_end.tv_sec - t_ddMC_start.tv_sec + 1e-6*(t_ddMC_end.tv_usec - t_ddMC_start.tv_usec));
    1265 #endif
     1253#endif 
    12661254                /* We loop through each row of P normalize it by making all
    12671255                * entries integer ones and add the resulting vector to the
    12681256                * int matrix facet::codim2Facets */
    12691257                for (int jj=1;jj<=/*ddakt*/P->rowsize;jj++)
    1270                 {
     1258                {                                       
    12711259                        fAct->numCodim2Facets++;
    1272                         if(fAct->numCodim2Facets==1)
    1273                         {
    1274                                 fAct->codim2Ptr = new facet(2);
     1260                        if(fAct->numCodim2Facets==1)                                   
     1261                        {                                               
     1262                                fAct->codim2Ptr = new facet(2);                                         
    12751263                                codim2Act = fAct->codim2Ptr;
    12761264                        }
     
    12971285#endif
    12981286                        codim2Act->setFacetNormal(n);
    1299                         delete n;
    1300                 }
     1287                        delete n;                                       
     1288                }               
    13011289                /*We check whether the facet spanned by the codim-2 facets
    13021290                * intersects with the positive orthant. Otherwise we define this
    1303                 * facet to be non-flippable. Works since we set the appropriate
     1291                * facet to be non-flippable. Works since we set the appropriate 
    13041292                * linearity for ddakt above.
    13051293                */
    1306                 TODO It might be faster to compute jus the implied equations instead of a relative interior point
    1307                 int64vec *iv_intPoint = new int64vec(this->numVars);
    1308                 dd_MatrixPtr shiftMatrix;
    1309                 dd_MatrixPtr intPointMatrix;
    1310                 shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
    1311                 for(int kk=0;kk<this->numVars;kk++)
    1312                 {
    1313                         dd_set_si(shiftMatrix->matrix[kk][0],1);
    1314                         dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
    1315                 }
    1316                 intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
    1317  #ifdef gfanp
    1318                 timeval t_iP_start, t_iP_end;
    1319                 gettimeofday(&t_iP_start, 0);
    1320  #endif
    1321                 interiorPoint(intPointMatrix,*iv_intPoint);
    1322  //             dd_rowset impl_linste,lbasis;
    1323  //             dd_LPSolutionPtr lps=NULL;
    1324  //             dd_ErrorType err;
    1325  //             dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
    1326  #ifdef gfanp
    1327                 gettimeofday(&t_iP_end, 0);
    1328                 t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
    1329  #endif
    1330                 for(int ll=0;ll<this->numVars;ll++)
    1331                 {
    1332                         if( (*iv_intPoint)[ll] < 0 )
    1333                         {
    1334                                 fAct->isFlippable=FALSE;
    1335                                 break;
    1336                         }
    1337                 }
     1294                //TODO It might be faster to compute jus the implied equations instead of a relative interior point
     1295//              int64vec *iv_intPoint = new int64vec(this->numVars);
     1296//              dd_MatrixPtr shiftMatrix;
     1297//              dd_MatrixPtr intPointMatrix;
     1298//              shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
     1299//              for(int kk=0;kk<this->numVars;kk++)
     1300//              {
     1301//                      dd_set_si(shiftMatrix->matrix[kk][0],1);
     1302//                      dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
     1303//              }
     1304//              intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
     1305// #ifdef gfanp
     1306//              timeval t_iP_start, t_iP_end;
     1307//              gettimeofday(&t_iP_start, 0);
     1308// #endif               
     1309//              interiorPoint(intPointMatrix,*iv_intPoint);
     1310// //           dd_rowset impl_linste,lbasis;
     1311// //           dd_LPSolutionPtr lps=NULL;
     1312// //           dd_ErrorType err;
     1313// //           dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
     1314// #ifdef gfanp
     1315//              gettimeofday(&t_iP_end, 0);
     1316//              t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
     1317// #endif
     1318//              for(int ll=0;ll<this->numVars;ll++)
     1319//              {
     1320//                      if( (*iv_intPoint)[ll] < 0 )
     1321//                      {
     1322//                              fAct->isFlippable=FALSE;
     1323//                              break;
     1324//                      }
     1325//              }
    13381326                /*End of check*/
    13391327                /*This test should be way less time consuming*/
     
    13581346                }
    13591347                if(containsStrictlyPosRay==FALSE)
    1360                 TODO Not sufficient. Intersect with pos orthant for pos int
     1348                //TODO Not sufficient. Intersect with pos orthant for pos int
    13611349                        fAct->isFlippable=FALSE;
    13621350#ifdef gfanp
     
    13651353#endif
    13661354                /**/
    1367                 fAct = fAct->next;
     1355                fAct = fAct->next;     
    13681356                dd_FreeMatrix(ddakt);
    13691357                dd_FreeMatrix(P);
    1370         }for
     1358        }//for 
    13711359        dd_FreeMatrix(ddineq);
    13721360#ifdef gfanp
     
    13751363#endif
    13761364}
    1377 
     1365               
    13781366/** Really extremal rays this time ;)
    13791367* Extremal rays are unique modulo the homogeneity space.
     
    13831371* checking whether the inner product of the ray with the normal is zero.
    13841372* We use ivAdd here which returns a new int64vec. Therefore we need to avoid
    1385 * a memory leak which would be cause by the line
     1373* a memory leak which would be cause by the line 
    13861374* iv=ivAdd(iv,b)
    1387 * So we keep pointer tmp to iv and delete(tmp), so there should not occur a
     1375* So we keep pointer tmp to iv and delete(tmp), so there should not occur a 
    13881376* memleak
    13891377* TODO normalization
     
    13971385        gettimeofday(&poly_start,0);
    13981386#endif
    1399         Add lineality space - dd_LinealitySpace
     1387        //Add lineality space - dd_LinealitySpace
    14001388        dd_MatrixPtr ddineq;
    1401         dd_ErrorType err;
    1402         if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
    1403                 ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
    1404         else
    1405                 ddineq = dd_CopyMatrix(gc.ddFacets);
     1389        dd_ErrorType err;       
     1390//      if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
     1391//              ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
     1392//      else
     1393//              ddineq = dd_CopyMatrix(gc.ddFacets);
    14061394        ddineq = (dd_LinealitySpace->rowsize>0) ? dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace) : dd_CopyMatrix(gc.ddFacets);
    14071395        /* In case the input is non-homogeneous we add constrains for the positive orthant.
    1408         * This is justified by the fact that for non-homog ideals we only consider the
     1396        * This is justified by the fact that for non-homog ideals we only consider the 
    14091397        * restricted fan. This way we can be sure to find strictly positive interior points.
    14101398        * This in turn makes life easy when checking for flippability!
     
    14211409                assert(ddineq);
    14221410                dd_FreeMatrix(ddPosRestr);
    1423         }
     1411        }       
    14241412        dd_PolyhedraPtr ddPolyh;
    14251413        ddPolyh = dd_DDMatrix2Poly(ddineq, &err);
     
    14341422        /* Compute interior point on the fly*/
    14351423        int64vec *ivIntPointOfCone = new int64vec(this->numVars);
    1436         mpq_t *colSum = new mpq_t[this->numVars];
    1437         int denom[this->numVars];//denominators of colSum
    1438         NOTE TODO need to gcd of rows and factor out! -> makeInt
    1439         /*for(int jj=0;jj<this->numVars;jj++)
    1440         {
    1441                 mpq_init(colSum[jj]);
    1442                 for(int ii=0;ii<P->rowsize;ii++)
    1443                 {
    1444                         mpq_t tmp; mpq_init(tmp);
    1445                         mpq_t sum; mpq_init(sum);
    1446                         mpq_set(sum,colSum[jj]);
    1447                         mpq_add(tmp,sum,P->matrix[ii][jj+1]);
    1448                         mpq_set(colSum[jj],tmp);
    1449                         mpq_clear(tmp);
    1450                         mpq_clear(sum);
    1451                 }
    1452                 mpz_t den; mpz_init(den);
    1453                 mpq_get_den(den,colSum[jj]);
    1454                 denom[jj]=(int)mpz_get_si(den);
    1455                 mpz_clear(den);
    1456         }
    1457         Now compute lcm of denominators of colSum[jj]
    1458         NOTE gcd as well and normalise instantly?
    1459         mpz_t kgV; mpz_init(kgV);
    1460         mpz_set_str(kgV,"1",10);
    1461         mpz_t den; mpz_init(den);
    1462         mpz_t tmp; mpz_init(tmp);
    1463         mpq_get_den(tmp,colSum[0]);
    1464         for (int ii=0;ii<(this->numVars)-1;ii++)
    1465         {
    1466                 mpq_get_den(den,colSum[ii+1]);
    1467                 mpz_lcm(kgV,tmp,den);
    1468                 mpz_set(tmp, kgV);
    1469         }
    1470         mpq_t qkgV;
    1471         mpq_init(qkgV);
    1472         mpq_set_z(qkgV,kgV);
    1473         mpz_clear(kgV);
    1474         mpz_clear(den);
    1475         mpz_clear(tmp);*/
    14761424        int64vec *foo = new int64vec(this->numVars);
    14771425        for(int ii=0;ii<P->rowsize;ii++)
    14781426        {
    1479                 int64vec *foo = new int64vec(this->numVars);
     1427//              int64vec *foo = new int64vec(this->numVars);
    14801428                int64vec *tmp = ivIntPointOfCone;
    14811429                makeInt(P,ii+1,*foo);
    14821430                ivIntPointOfCone = iv64Add(ivIntPointOfCone,foo);
    14831431                delete tmp;
    1484                 delete foo;
     1432//              delete foo;
    14851433        }
    14861434        delete foo;
     
    14881436        for (int ii=0;ii<(this->numVars);ii++)
    14891437        {
    1490                 mpq_t product;
    1491                 mpq_init(product);
    1492                 mpq_mul(product,qkgV,colSum[ii]);
    1493                 (*ivIntPointOfCone)[ii]=(int64)mpz_get_d(mpq_numref(product));
    1494                 if( (*ivIntPointOfCone)[ii]>INT_MAX )
     1438//              mpq_t product;
     1439//              mpq_init(product);
     1440//              mpq_mul(product,qkgV,colSum[ii]);
     1441//              (*ivIntPointOfCone)[ii]=(int64)mpz_get_d(mpq_numref(product));
     1442                if( (*ivIntPointOfCone)[ii]>INT_MAX ) 
    14951443                        WarnS("Interior point exceeds INT_MAX!\n");
    1496                 mpq_clear(product);
    1497                 Compute intgcd
     1444//              mpq_clear(product);
     1445                //Compute intgcd
    14981446                ggT=int64gcd(ggT,(*ivIntPointOfCone)[ii]);
    14991447        }
    1500 
    1501         Divide out a common gcd > 1
     1448       
     1449        //Divide out a common gcd > 1
    15021450        if(ggT>1)
    15031451        {
     
    15081456                }
    15091457        }
    1510         mpq_clear(qkgV);
    1511         delete [] colSum;
     1458//      mpq_clear(qkgV);
     1459//      delete [] colSum;
    15121460        /*For homogeneous input (like Det3,3,5) the int points may be negative. So add a suitable multiple of (1,_,1)*/
    15131461        if(hasHomInput==TRUE && iv64isStrictlyPositive(ivIntPointOfCone)==FALSE)
     
    15171465                for(int ii=0;ii<this->numVars;ii++)
    15181466                {
    1519                         (*ivOne)[ii]=1;
     1467//                      (*ivOne)[ii]=1;
    15201468                        if ((*ivIntPointOfCone)[ii]<maxNegEntry) maxNegEntry=(*ivIntPointOfCone)[ii];
    15211469                }
    15221470                maxNegEntry *= -1;
    1523                 maxNegEntry++;To be on the safe side
     1471                maxNegEntry++;//To be on the safe side
    15241472                for(int ii=0;ii<this->numVars;ii++)
    15251473                        (*ivOne)[ii]=maxNegEntry;
     
    15271475                ivIntPointOfCone=iv64Add(ivIntPointOfCone,ivOne);
    15281476                delete(tmp);
    1529                 while( !iv64isStrictlyPositive(ivIntPointOfCone) )
    1530                 {
    1531                         int64vec *tmp = ivIntPointOfCone;
    1532                         for(int jj=0;jj<this->numVars;jj++)
    1533                                 (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
    1534                         ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
    1535                         delete tmp;
    1536                 }
     1477//              while( !iv64isStrictlyPositive(ivIntPointOfCone) )
     1478//              {
     1479//                      int64vec *tmp = ivIntPointOfCone;
     1480//                      for(int jj=0;jj<this->numVars;jj++)
     1481//                              (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
     1482//                      ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
     1483//                      delete tmp;                             
     1484//              }
    15371485                delete ivOne;
    15381486                int64 ggT=(*ivIntPointOfCone)[0];
     
    15451493                }
    15461494        }
    1547         assert(iv64isStrictlyPositive(ivIntPointOfCone));
    1548 
     1495//      assert(iv64isStrictlyPositive(ivIntPointOfCone));
     1496       
    15491497        this->setIntPoint(ivIntPointOfCone);
    15501498        delete(ivIntPointOfCone);
    15511499        /* end of interior point computation*/
    1552 
    1553         Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
     1500       
     1501        //Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
    15541502        int rows=P->rowsize;
    15551503        facet *fAct=gc.facetPtr;
    1556         Construct an array to hold the extremal rays of the cone
    1557         this->gcRays = (int64vec**)omAlloc0(sizeof(int64vec*)*P->rowsize);
     1504        //Construct an array to hold the extremal rays of the cone
     1505        this->gcRays = (int64vec**)omAlloc0(sizeof(int64vec*)*P->rowsize);     
    15581506        for(int ii=0;ii<P->rowsize;ii++)
    15591507        {
    15601508                int64vec *rowvec = new int64vec(this->numVars);
    1561                 makeInt(P,ii+1,*rowvec);get an integer entry instead of rational, rowvec is primitve
     1509                makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve
    15621510                this->gcRays[ii] = iv64Copy(rowvec);
    15631511                delete rowvec;
    1564         }
     1512        }       
    15651513        this->numRays=P->rowsize;
    1566         Check which rays belong to which facet
     1514        //Check which rays belong to which facet
    15671515        while(fAct!=NULL)
    15681516        {
    1569                 const int64vec *fNormal; = new int64vec(this->numVars);
    1570                 fNormal = fAct->getRef2FacetNormal();->getFacetNormal();
     1517                const int64vec *fNormal;// = new int64vec(this->numVars);
     1518                fNormal = fAct->getRef2FacetNormal();//->getFacetNormal();
    15711519                int64vec *ivIntPointOfFacet = new int64vec(this->numVars);
    15721520                for(int ii=0;ii<rows;ii++)
    1573                 {
     1521                {                       
    15741522                        if(dotProduct(*fNormal,this->gcRays[ii])==0)
    15751523                        {
    1576                                 int64vec *tmp = ivIntPointOfFacet;Prevent memleak
     1524                                int64vec *tmp = ivIntPointOfFacet;//Prevent memleak
    15771525                                fAct->numCodim2Facets++;
    15781526                                facet *codim2Act;
    1579                                 if(fAct->numCodim2Facets==1)
    1580                                 {
    1581                                         fAct->codim2Ptr = new facet(2);
     1527                                if(fAct->numCodim2Facets==1)                                   
     1528                                {                                               
     1529                                        fAct->codim2Ptr = new facet(2);                                         
    15821530                                        codim2Act = fAct->codim2Ptr;
    15831531                                }
     
    15871535                                        codim2Act = codim2Act->next;
    15881536                                }
    1589                                 codim2Act->setFacetNormal(rowvec);
    1590                                 Rather just let codim2Act point to the corresponding int64vec of gcRays
     1537                                //codim2Act->setFacetNormal(rowvec);
     1538                                //Rather just let codim2Act point to the corresponding int64vec of gcRays
    15911539                                codim2Act->fNormal=this->gcRays[ii];
    1592                                 fAct->numRays++;
    1593                                 Memleak avoided via tmp
     1540                                fAct->numRays++;                                 
     1541                                //Memleak avoided via tmp
    15941542                                ivIntPointOfFacet=iv64Add(ivIntPointOfFacet,this->gcRays[ii]);
    1595                                 Now tmp still points to the OLD address of ivIntPointOfFacet
     1543                                //Now tmp still points to the OLD address of ivIntPointOfFacet
    15961544                                delete(tmp);
    1597 
     1545                                       
    15981546                        }
    1599                 }For non-homog input ivIntPointOfFacet should already be >0 here
    1600                 if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
    1601                 if we have no strictly pos ray but the input is homogeneous
    1602                 then add a suitable multiple of (1,...,1)
     1547                }//For non-homog input ivIntPointOfFacet should already be >0 here
     1548//              if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
     1549                //if we have no strictly pos ray but the input is homogeneous
     1550                //then add a suitable multiple of (1,...,1)
    16031551                if( !iv64isStrictlyPositive(ivIntPointOfFacet) && hasHomInput==TRUE)
    16041552                {
     
    16111559                                for(int jj=0;jj<this->numVars;jj++)
    16121560                                {
    1613                                         (*ivOne)[jj] = (*ivOne)[jj] << 1; times 2
     1561                                        (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
    16141562                                }
    16151563                                ivIntPointOfFacet = iv64Add(ivIntPointOfFacet/*diff*/,ivOne);
    1616                                 delete tmp;
     1564                                delete tmp;                             
    16171565                        }
    16181566                        delete ivOne;
     
    16251573                        for(int ii=0;ii<this->numVars;ii++)
    16261574                                (*ivIntPointOfFacet)[ii] /= ggT;
    1627                 }
     1575                }                       
    16281576                fAct->setInteriorPoint(ivIntPointOfFacet);
    1629 
     1577               
    16301578                delete(ivIntPointOfFacet);
    1631                 Now (if we have at least 3 variables) do a bubblesort on the rays
     1579                //Now (if we have at least 3 variables) do a bubblesort on the rays
    16321580                /*if(this->numVars>2)
    16331581                {
     
    16431591                        do
    16441592                        {
    1645                                 exchanged=FALSE;n=fAct->numRays-1;
     1593                                exchanged=FALSE;//n=fAct->numRays-1;                           
    16461594                                for(unsigned ii=0;ii<=n-1;ii++)
    1647                                 {
     1595                                {                                       
    16481596                                        if((A[ii]->fNormal)->compare((A[ii+1]->fNormal))==1)
    16491597                                        {
    1650                                                 Swap rays
     1598                                                //Swap rays
    16511599                                                cout << "Swapping ";
    16521600                                                A[ii]->fNormal->show(1,0); cout << " with "; A[ii+1]->fNormal->show(1,0); cout << endl;
     
    16571605                                                if(ii==0)
    16581606                                                        fAct->codim2Ptr=A[ii+1];
    1659                                                 end swap
    1660                                                 facet *tmp=A[ii];swap in list
     1607                                                //end swap
     1608                                                facet *tmp=A[ii];//swap in list
    16611609                                                A[ii+1]=A[ii];
    16621610                                                A[ii]=tmp;
    1663                                                 tmp=NULL;
    1664                                         }
     1611//                                              tmp=NULL;                       
     1612                                        }                                       
    16651613                                }
    1666                                 n--;
     1614                                n--;                   
    16671615                        }while(exchanged==TRUE && n>=0);
    1668                 }*/if pVariables>2
    1669                 delete fNormal;
     1616                }*///if pVariables>2
     1617//              delete fNormal;         
    16701618                fAct = fAct->next;
    1671         }end of facet checking
     1619        }//end of facet checking
    16721620        dd_FreeMatrix(P);
    1673         Now all extremal rays should be set w.r.t their respective fNormal
    1674         TODO Not sufficient -> vol2 II/125&127
    1675         NOTE Sufficient according to cddlibs doc. These ARE rays
     1621        //Now all extremal rays should be set w.r.t their respective fNormal
     1622        //TODO Not sufficient -> vol2 II/125&127
     1623        //NOTE Sufficient according to cddlibs doc. These ARE rays
    16761624        //What the hell... let's just take interior points
    16771625        if(gcone::hasHomInput==FALSE)
     
    16801628                while(fAct!=NULL)
    16811629                {
    1682                         bool containsStrictlyPosRay=FALSE;
    1683                         facet *codim2Act;
    1684                         codim2Act = fAct->codim2Ptr;
    1685                         while(codim2Act!=NULL)
    1686                         {
    1687                                 int64vec *rayvec;
    1688                                 rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
    1689                                 //int negCtr=0;
    1690                                 if(iv64isStrictlyPositive(rayvec))
    1691                                 {
    1692                                         containsStrictlyPosRay=TRUE;
    1693                                         delete(rayvec);
    1694                                         break;
    1695                                 }
    1696                                 delete(rayvec);
    1697                                 codim2Act = codim2Act->next;
    1698                         }
    1699                         if(containsStrictlyPosRay==FALSE)
    1700                                 fAct->isFlippable=FALSE;
     1630//                      bool containsStrictlyPosRay=FALSE;
     1631//                      facet *codim2Act;
     1632//                      codim2Act = fAct->codim2Ptr;
     1633//                      while(codim2Act!=NULL)
     1634//                      {
     1635//                              int64vec *rayvec;
     1636//                              rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
     1637//                              //int negCtr=0;
     1638//                              if(iv64isStrictlyPositive(rayvec))
     1639//                              {
     1640//                                      containsStrictlyPosRay=TRUE;
     1641//                                      delete(rayvec);
     1642//                                      break;
     1643//                              }                               
     1644//                              delete(rayvec);
     1645//                              codim2Act = codim2Act->next;
     1646//                      }
     1647//                      if(containsStrictlyPosRay==FALSE)
     1648//                              fAct->isFlippable=FALSE;
    17011649                        if(!iv64isStrictlyPositive(fAct->interiorPoint))
    17021650                                fAct->isFlippable=FALSE;
    17031651                        fAct = fAct->next;
    17041652                }
    1705         }hasHomInput?
     1653        }//hasHomInput?
    17061654#ifdef gfanp
    17071655        gettimeofday(&end, 0);
    17081656        t_getExtremalRays += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    1709 #endif
     1657#endif 
    17101658}
    17111659
     
    17131661void gcone::orderRays()
    17141662{
    1715    qsort(gcRays,sizeof(int64vec),int64vec::compare);
    1716 }
    1717 
     1663//   qsort(gcRays,sizeof(int64vec),int64vec::compare);
     1664}
     1665       
    17181666inline bool gcone::iv64isStrictlyPositive(const int64vec * iv64)
    17191667{
     
    17251673                        res=FALSE;
    17261674                        break;
    1727                 }
     1675                }               
    17281676        }
    17291677        return res;
    17301678}
    17311679
    1732 /** \brief Compute the Groebner Basis on the other side of a shared facet
     1680/** \brief Compute the Groebner Basis on the other side of a shared facet 
    17331681 *
    17341682 * Implements algorithm 4.3.2 from Anders' thesis.
     
    17381686 * Parallelity is checked using basic linear algebra. See gcone::isParallel.
    17391687 * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
    1740  * computing an interior point of the facet and taking all terms having the same weight with respect
     1688 * computing an interior point of the facet and taking all terms having the same weight with respect 
    17411689 * to this interior point.
    17421690 *\param ideal, facet
    17431691 * Input: a marked,reduced Groebner basis and a facet
    17441692 */
    1745 inline void gcone::flip(ideal gb, facet *f)             Compute "the other side"
    1746 {
     1693inline void gcone::flip(ideal gb, facet *f)             //Compute "the other side"
     1694{       
    17471695#ifdef gfanp
    17481696        timeval start, end;
    17491697        gettimeofday(&start, 0);
    1750 #endif
    1751         int64vec *fNormal; = new int64vec(this->numVars);       //facet normal, check for parallelity
    1752         fNormal = f->getFacetNormal();  read this->fNormal;
     1698#endif         
     1699        int64vec *fNormal;// = new int64vec(this->numVars);     //facet normal, check for parallelity                   
     1700        fNormal = f->getFacetNormal();  //read this->fNormal;
    17531701#ifdef gfan_DEBUG
    1754         std::cout << "running gcone::flip" << std::endl;
     1702//      std::cout << "running gcone::flip" << std::endl;
    17551703        printf("flipping UCN %i over facet",this->getUCN());
    17561704        fNormal->show(1,0);
    1757         printf(") with UCN %i\n",f->getUCN() );
     1705        printf(") with UCN %i\n",f->getUCN() ); 
    17581706#endif
    17591707        if(this->getUCN() != f->getUCN())
     
    17631711        }
    17641712        /*1st step: Compute the initial ideal*/
    1765         /*poly initialFormElement[IDELEMS(gb)];*/       array of #polys in GB to store initial form
     1713        /*poly initialFormElement[IDELEMS(gb)];*/       //array of #polys in GB to store initial form
    17661714        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
    1767 
     1715       
    17681716        computeInv(gb,initialForm,*fNormal);
    17691717
     
    17711719/*      cout << "Initial ideal is: " << endl;
    17721720        idShow(initialForm);
    1773         f->printFlipGB();*/
    1774         cout << "===" << endl;
    1775 #endif
     1721        //f->printFlipGB();*/
     1722//      cout << "===" << endl;
     1723#endif                 
    17761724        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
    17771725        /*Substep 2.1
    1778         compute $G_{-\alpha}(in_v(I))
     1726        compute $G_{-\alpha}(in_v(I)) 
    17791727        see journal p. 66
    17801728        NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
     
    17831731        ring srcRing=currRing;
    17841732        ring tmpRing;
    1785 
     1733                       
    17861734        if( (srcRing->order[0]!=ringorder_a))
    17871735        {
    1788                 int64vec *iv; = new int64vec(this->numVars);
    1789                 iv = ivNeg(fNormal);ivNeg uses iv64Copy -> new
    1790                 tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
     1736                int64vec *iv;// = new int64vec(this->numVars);
     1737                iv = ivNeg(fNormal);//ivNeg uses iv64Copy -> new
     1738//              tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
    17911739                tmpRing=rCopyAndAddWeight(srcRing,iv);
    17921740                delete iv;
     
    18051753                tmpRing->block1[0]=length;
    18061754                rComplete(tmpRing);
    1807                 omFree(A);
    1808         }
    1809         delete fNormal;
    1810         rChangeCurrRing(tmpRing);
    1811 
    1812         ideal ina;
     1755                //omFree(A);
     1756        }
     1757        delete fNormal; 
     1758        rChangeCurrRing(tmpRing);       
     1759                       
     1760        ideal ina;                     
    18131761        ina=idrCopyR(initialForm,srcRing);
    18141762        idDelete(&initialForm);
    18151763        ideal H;
    1816         H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
     1764//      H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
    18171765#ifdef gfanp
    18181766        timeval t_kStd_start, t_kStd_end;
     
    18221770                H=kStd(ina,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
    18231771        else
    1824                 H=kStd(ina,NULL,isNotHomog,NULL);       This is \mathcal(G)_{>_-\alpha}(in_v(I))
     1772                H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
    18251773#ifdef gfanp
    18261774        gettimeofday(&t_kStd_end, 0);
     
    18351783        rChangeCurrRing(srcRing);
    18361784        ideal srcRing_H;
    1837         ideal srcRing_HH;
     1785        ideal srcRing_HH;                       
    18381786        srcRing_H=idrCopyR(H,tmpRing);
    18391787        //H is needed further below, so don't idDelete here
    18401788        srcRing_HH=ffG(srcRing_H,this->gcBasis);
    18411789        idDelete(&srcRing_H);
    1842 
     1790               
    18431791        /*Substep 2.2.1
    18441792         * Mark according to G_-\alpha
     
    18461794         * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
    18471795         * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
    1848          * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
     1796         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we 
    18491797         * compute the difference accordingly
    18501798        */
     
    18521800        timeval t_markings_start, t_markings_end;
    18531801        gettimeofday(&t_markings_start, 0);
    1854 #endif
     1802#endif         
    18551803        bool markingsAreCorrect=FALSE;
    18561804        dd_MatrixPtr intPointMatrix;
    18571805        int iPMatrixRows=0;
    1858         dd_rowrange aktrow=0;
     1806        dd_rowrange aktrow=0;                   
    18591807        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
    18601808        {
    18611809                poly aktpoly=(poly)srcRing_HH->m[ii];//This is a pointer, so don't pDelete
    1862                 iPMatrixRows = iPMatrixRows+pLength(aktpoly);
     1810                iPMatrixRows = iPMatrixRows+pLength(aktpoly);           
    18631811        }
    18641812        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
    18651813         * construction of the differences
    18661814         */
    1867         intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
    1868         intPointMatrix->numbtype=dd_Integer;    NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
    1869 
     1815        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); 
     1816        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
     1817       
    18701818        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
    18711819        {
    1872                 markingsAreCorrect=FALSE;       crucial to initialise here
     1820                markingsAreCorrect=FALSE;       //crucial to initialise here
    18731821                poly aktpoly=srcRing_HH->m[ii]; //Only a pointer, so don't pDelete
    18741822                /*Comparison of leading monomials is done via exponent vectors*/
     
    18781826                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
    18791827                        pGetExpV(aktpoly,src_ExpV);
    1880                         rChangeCurrRing(tmpRing);       this ring change is crucial!
     1828                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
    18811829                        poly p=pCopy(H->m[ii]);
    18821830                        pGetExpV(p/*pCopy(H->m[ii])*/,dst_ExpV);
     
    18871835                        {
    18881836#ifdef gfan_DEBUG
    1889                                 cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
     1837//                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
    18901838#endif
    18911839                                if (src_ExpV[kk]!=dst_ExpV[kk])
     
    18931841                                        expVAreEqual=FALSE;
    18941842                                }
    1895                         }
     1843                        }                                                               
    18961844                        if (expVAreEqual==TRUE)
    18971845                        {
    1898                                 markingsAreCorrect=TRUE; everything is fine
     1846                                markingsAreCorrect=TRUE; //everything is fine
    18991847#ifdef gfan_DEBUG
    1900                                 cout << "correct markings" << endl;
    1901 #endif
    1902                         }if (pHead(aktpoly)==pHead(H->m[jj])
     1848//                              cout << "correct markings" << endl;
     1849#endif
     1850                        }//if (pHead(aktpoly)==pHead(H->m[jj])
    19031851                        omFree(src_ExpV);
    19041852                        omFree(dst_ExpV);
    1905                 }for (int jj=0;jj<IDELEMS(H);jj++)
    1906 
     1853                }//for (int jj=0;jj<IDELEMS(H);jj++)
     1854               
    19071855                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
    19081856                if (markingsAreCorrect==TRUE)
     
    19131861                {
    19141862                        rChangeCurrRing(tmpRing);
    1915                         pGetExpV(pHead(H->m[ii]),leadExpV); We use H->m[ii] as leading monomial
     1863                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
    19161864                        rChangeCurrRing(srcRing);
    19171865                }
    1918                 /*compute differences of the expvects*/
     1866                /*compute differences of the expvects*/                         
    19191867                while (pNext(aktpoly)!=NULL)
    19201868                {
    19211869                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    1922                         /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
     1870                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term) 
    19231871                        is not omitted when computing the differences*/
    19241872                        if(markingsAreCorrect==TRUE)
     
    19341882                        int ctr=0;
    19351883                        for (int jj=0;jj<this->numVars;jj++)
    1936                         {
    1937                                 /*Store into ddMatrix*/
     1884                        {                               
     1885                                /*Store into ddMatrix*/                         
    19381886                                if(leadExpV[jj+1]-v[jj+1])
    19391887                                        ctr++;
     
    19411889                        }
    19421890                        /*It ought to be more sensible to avoid 0-rows in the first place*/
    1943                         if(ctr==this->numVars)//We have a 0-row
    1944                                 dd_MatrixRowRemove(&intPointMatrix,aktrow);
    1945                         else
     1891//                      if(ctr==this->numVars)//We have a 0-row
     1892//                              dd_MatrixRowRemove(&intPointMatrix,aktrow);
     1893//                      else
    19461894                                aktrow +=1;
    19471895                        omFree(v);
    19481896                }
    19491897                omFree(leadExpV);
    1950         }for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
     1898        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
    19511899#ifdef gfanp
    19521900        gettimeofday(&t_markings_end, 0);
     
    19601908        preprocessInequalities(intPointMatrix);
    19611909        /*Now we add the constraint for the standard simplex*/
    1962         dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
    1963         for (int jj=1;jj<=this->numVars;jj++)
    1964         {
    1965                 dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
    1966         }
     1910//      dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
     1911//      for (int jj=1;jj<=this->numVars;jj++)
     1912//      {
     1913//              dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
     1914//      }
    19671915        //Let's make sure we compute interior points from the positive orthant
    1968         dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
    1969 
    1970         int jj=1;
    1971         for (int ii=0;ii<this->numVars;ii++)
    1972         {
    1973                 dd_set_si(posRestr->matrix[ii][jj],1);
    1974                 jj++;
    1975         }
     1916//      dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
     1917//     
     1918//      int jj=1;
     1919//      for (int ii=0;ii<this->numVars;ii++)
     1920//      {
     1921//              dd_set_si(posRestr->matrix[ii][jj],1);
     1922//              jj++;                                                   
     1923//      }
    19761924        /*We create a matrix containing the standard simplex
    19771925        * and constraints to assure a strictly positive point
     
    19831931                {
    19841932                        dd_set_si(posRestr->matrix[ii][0],-1);
    1985                         for(int jj=1;jj<=this->numVars;jj++)
    1986                                 dd_set_si(posRestr->matrix[ii][jj],1);
     1933                        for(int jj=1;jj<=this->numVars;jj++)                   
     1934                                dd_set_si(posRestr->matrix[ii][jj],1);                 
    19871935                }
    19881936                else
     
    20051953        dd_rowindex newpos;
    20061954
    2007         NOTE Here we should remove interiorPoint and instead
    2008         create and ordering like (a(omega),a(fNormal),dp)
    2009         if(this->ivIntPt==NULL)
    2010                 interiorPoint(intPointMatrix, *iv_weight);      iv_weight now contains the interior point
    2011         else
    2012                 iv_weight=this->getIntPoint();
     1955        //NOTE Here we should remove interiorPoint and instead
     1956        // create and ordering like (a(omega),a(fNormal),dp)
     1957//      if(this->ivIntPt==NULL)
     1958                interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
     1959//      else
     1960//              iv_weight=this->getIntPoint();
    20131961        dd_FreeMatrix(intPointMatrix);
    20141962        /*Crude attempt for interior point */
     
    20341982        gettimeofday(&t_dd_end, 0);
    20351983        t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec));
    2036 #endif
    2037 
     1984#endif 
     1985       
    20381986        /*Step 3
    2039         * turn the minimal basis into a reduced one */
    2040         NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
    2041          Thus:
    2042         ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
     1987        * turn the minimal basis into a reduced one */                 
     1988        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
     1989        // Thus:
     1990        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
    20431991        ring dstRing=rCopy0(tmpRing);
    20441992        int length=iv_weight->length();
     
    20542002        delete iv_weight;
    20552003
    2056         ideal dstRing_I;
     2004        ideal dstRing_I;                       
    20572005        dstRing_I=idrCopyR(srcRing_HH,srcRing);
    2058         idDelete(&srcRing_HH); Hmm.... causes trouble - no more
    2059         dstRing_I=idrCopyR(inputIdeal,srcRing);
     2006        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more
     2007        //dstRing_I=idrCopyR(inputIdeal,srcRing);
    20602008        BITSET save=test;
    20612009        test|=Sy_bit(OPT_REDSB);
    20622010        test|=Sy_bit(OPT_REDTAIL);
    20632011#ifdef gfan_DEBUG
    2064         test|=Sy_bit(6);        //OPT_DEBUG
     2012//      test|=Sy_bit(6);        //OPT_DEBUG
    20652013#endif
    20662014        ideal tmpI;
    2067         NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
    2068         tmpI = idrCopyR(this->inputIdeal,this->baseRing);
     2015        //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
     2016        //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
    20692017        tmpI = dstRing_I;
    20702018#ifdef gfanp
     
    20802028#endif
    20812029        idDelete(&tmpI);
    2082         idNorm(dstRing_I);
    2083         kInterRed(dstRing_I);
     2030        idNorm(dstRing_I);                     
     2031//      kInterRed(dstRing_I);
    20842032        idSkipZeroes(dstRing_I);
    20852033        test=save;
    20862034        /*End of step 3 - reduction*/
    2087 
    2088         f->setFlipGB(dstRing_I);store the flipped GB
    2089         idDelete(&dstRing_I);
    2090         f->flipRing=rCopy(dstRing);     store the ring on the other side
     2035                       
     2036        f->setFlipGB(dstRing_I);//store the flipped GB
     2037//      idDelete(&dstRing_I);
     2038        f->flipRing=rCopy(dstRing);     //store the ring on the other side
    20912039#ifdef gfan_DEBUG
    20922040        printf("Flipped GB is UCN %i:\n",counter+1);
    20932041        idDebugPrint(dstRing_I);
    20942042        printf("\n");
    2095 #endif
    2096         idDelete(&dstRing_I);
    2097         rChangeCurrRing(srcRing);       return to the ring we started the computation of flipGB in
     2043#endif 
     2044        idDelete(&dstRing_I);   
     2045        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
    20982046        rDelete(dstRing);
    20992047#ifdef gfanp
     
    21012049        time_flip += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    21022050#endif
    2103 }void flip(ideal gb, facet *f)
     2051}//void flip(ideal gb, facet *f)
    21042052
    21052053/** \brief A slightly different approach to flipping
     
    21112059* will be from a ring with (a(),dp,C) and our resulting cone from (a(),a(),dp,C). Hence a way
    21122060* must be found to circumvent the sequence of a()'s growing to a ridiculous size.
    2113 * Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we
    2114 * do have an interior point of the cone by adding the extremal rays. So we replace
    2115 * the latter cone by a cone with (a(sum_of_rays),dp,C).
     2061* Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we 
     2062* do have an interior point of the cone by adding the extremal rays. So we replace 
     2063* the latter cone by a cone with (a(sum_of_rays),dp,C). 
    21162064* Con: It's incredibly ugly
    21172065* Pro: No messing around with readConeFromFile()
     
    21252073#endif
    21262074        const int64vec *fNormal;
    2127         fNormal = f->getRef2FacetNormal();/*->getFacetNormal();*/       read this->fNormal;
     2075        fNormal = f->getRef2FacetNormal();/*->getFacetNormal();*/       //read this->fNormal;
    21282076#ifdef gfan_DEBUG
    21292077        printf("flipping UCN %i over facet(",this->getUCN());
    21302078        fNormal->show(1,0);
    2131         printf(") with UCN %i\n",f->getUCN());
     2079        printf(") with UCN %i\n",f->getUCN()); 
    21322080#endif
    21332081        if(this->getUCN() != f->getUCN())
     
    21372085        }
    21382086        /*1st step: Compute the initial ideal*/
    2139         ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
     2087        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);     
    21402088        computeInv( gb, initialForm, *fNormal );
    21412089        ring srcRing=currRing;
    21422090        ring tmpRing;
    2143 
     2091       
    21442092        const int64vec *intPointOfFacet;
    21452093        intPointOfFacet=f->getInteriorPoint();
    2146         Now we need two blocks of ringorder_a!
    2147         May assume the same situation as in flip() here
     2094        //Now we need two blocks of ringorder_a!
     2095        //May assume the same situation as in flip() here                       
    21482096        if( (srcRing->order[0]!=ringorder_a/*64*/) && (srcRing->order[1]!=ringorder_a/*64*/) )
    21492097        {
    2150                 int64vec *iv = new int64vec(this->numVars);init with 1s, since we do not need a 2nd block here but later
    2151                 int64vec *iv_foo = new int64vec(this->numVars,1);//placeholder
    2152                 int64vec *ivw = ivNeg(const_cast<int64vec*>(fNormal));
     2098                int64vec *iv = new int64vec(this->numVars);//init with 1s, since we do not need a 2nd block here but later
     2099//              int64vec *iv_foo = new int64vec(this->numVars,1);//placeholder
     2100                int64vec *ivw = ivNeg(const_cast<int64vec*>(fNormal));         
    21532101                tmpRing=rCopyAndAddWeight2(srcRing,ivw/*intPointOfFacet*/,iv);
    21542102                delete iv;delete ivw;
    2155                 delete iv_foo;
    2156         }
    2157         else
     2103//              delete iv_foo;
     2104        }
     2105        else 
    21582106        {
    21592107                int64vec *iv=new int64vec(this->numVars);
     
    21682116                {
    21692117                        A1[jj] = -(*fNormal)[jj];
    2170                         A2[jj] = 1;-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
     2118                        A2[jj] = 1;//-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
    21712119                }
    21722120                omFree(tmpRing->wvhdl[0]);
    21732121                if(tmpRing->wvhdl[1]!=NULL)
    21742122                        omFree(tmpRing->wvhdl[1]);
    2175                 tmpRing->wvhdl[0]=(int*)A1;
     2123                tmpRing->wvhdl[0]=(int*)A1;             
    21762124                tmpRing->block1[0]=length;
    2177                 tmpRing->wvhdl[1]=(int*)A2;
     2125                tmpRing->wvhdl[1]=(int*)A2;             
    21782126                tmpRing->block1[1]=length;
    21792127                rComplete(tmpRing);*/
    21802128        }
    2181         delete fNormal; //NOTE Do not delete when using getRef2FacetNormal();
    2182         rChangeCurrRing(tmpRing);
    2183         Now currRing should have (a(),a(),dp,C)
    2184         ideal ina;
     2129//      delete fNormal; //NOTE Do not delete when using getRef2FacetNormal();
     2130        rChangeCurrRing(tmpRing);       
     2131        //Now currRing should have (a(),a(),dp,C)               
     2132        ideal ina;                     
    21852133        ina=idrCopyR(initialForm,srcRing);
    21862134        idDelete(&initialForm);
     
    21932141        test|=Sy_bit(OPT_REDSB);
    21942142        test|=Sy_bit(OPT_REDTAIL);
    2195         if(gcone::hasHomInput==TRUE)
     2143//      if(gcone::hasHomInput==TRUE)
    21962144                H=kStd(ina,NULL,testHomog/*isHomog*/,NULL/*,gcone::hilbertFunction*/);
    2197         else
    2198                 H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
     2145//      else
     2146//              H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
    21992147        test=save;
    22002148#ifdef gfanp
     
    22042152        idSkipZeroes(H);
    22052153        idDelete(&ina);
    2206 
     2154       
    22072155        rChangeCurrRing(srcRing);
    22082156        ideal srcRing_H;
    2209         ideal srcRing_HH;
     2157        ideal srcRing_HH;                       
    22102158        srcRing_H=idrCopyR(H,tmpRing);
    22112159        //H is needed further below, so don't idDelete here
    22122160        srcRing_HH=ffG(srcRing_H,this->gcBasis);
    22132161        idDelete(&srcRing_H);
    2214         Now BBA(srcRing_HH) with (a(),a(),dp)
     2162        //Now BBA(srcRing_HH) with (a(),a(),dp)
    22152163        /* Evil modification of currRing */
    22162164        ring dstRing=rCopy0(tmpRing);
     
    22222170        {
    22232171                A1[jj] = (*intPointOfFacet)[jj];
    2224                 A2[jj] = -(*ivw)[jj];TODO Or minus (*ivw)[ii] ??? NOTE minus
     2172                A2[jj] = -(*ivw)[jj];//TODO Or minus (*ivw)[ii] ??? NOTE minus
    22252173        }
    22262174        omFree(dstRing->wvhdl[0]);
    22272175        if(dstRing->wvhdl[1]!=NULL)
    22282176                omFree(dstRing->wvhdl[1]);
    2229         dstRing->wvhdl[0]=(int*)A1;
     2177        dstRing->wvhdl[0]=(int*)A1;             
    22302178        dstRing->block1[0]=length;
    2231         dstRing->wvhdl[1]=(int*)A2;
     2179        dstRing->wvhdl[1]=(int*)A2;             
    22322180        dstRing->block1[1]=length;
    22332181        rComplete(dstRing);
    22342182        rChangeCurrRing(dstRing);
    2235         ideal dstRing_I;
     2183        ideal dstRing_I;                       
    22362184        dstRing_I=idrCopyR(srcRing_HH,srcRing);
    2237         idDelete(&srcRing_HH); Hmm.... causes trouble - no more
     2185        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more       
    22382186        save=test;
    22392187        test|=Sy_bit(OPT_REDSB);
     
    22422190        tmpI = dstRing_I;
    22432191#ifdef gfanp
    2244         timeval t_kStd_start, t_kStd_end;
     2192//      timeval t_kStd_start, t_kStd_end;
    22452193        gettimeofday(&t_kStd_start,0);
    22462194#endif
    2247         if(gcone::hasHomInput==TRUE)
    2248                 dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
    2249         else
     2195//      if(gcone::hasHomInput==TRUE)
     2196//              dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
     2197//      else
    22502198                dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
    22512199#ifdef gfanp
     
    22542202#endif
    22552203        idDelete(&tmpI);
    2256         idNorm(dstRing_I);
     2204        idNorm(dstRing_I);                     
    22572205        idSkipZeroes(dstRing_I);
    22582206        test=save;
    22592207        /*End of step 3 - reduction*/
    2260 
     2208       
    22612209        f->setFlipGB(dstRing_I);
    22622210        f->flipRing=rCopy(dstRing);
    22632211        rDelete(tmpRing);
    22642212        rDelete(dstRing);
    2265         Now we should have dstRing with (a(),a(),dp,C)
    2266         This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
    2267         of cones in noRevS
     2213        //Now we should have dstRing with (a(),a(),dp,C)
     2214        //This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
     2215        //of cones in noRevS
    22682216        rChangeCurrRing(srcRing);
    22692217#ifdef gfanp
     
    22712219        time_flip2 += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    22722220#endif
    2273 }flip2
     2221}//flip2
    22742222
    22752223/** \brief Compute initial ideal
     
    22892237                poly aktpoly = (poly)gb->m[ii];//Ptr, so don't pDelete(aktpoly)
    22902238                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
    2291                 pGetExpV(aktpoly,leadExpV);     find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
     2239                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
    22922240                initialFormElement=pHead(aktpoly);
    2293                 int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
     2241//              int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    22942242                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
    22952243                {
    22962244                        int64vec *check = new int64vec(this->numVars);
    2297                         aktpoly=pNext(aktpoly); next term
     2245                        aktpoly=pNext(aktpoly); //next term
    22982246                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    2299                         pGetExpV(aktpoly,v);
     2247                        pGetExpV(aktpoly,v);           
    23002248                        /* Convert (int)v into (int64vec)check */
    2301                         bool notPar=FALSE;
     2249//                      bool notPar=FALSE;
    23022250                        for (int jj=0;jj<this->numVars;jj++)
    23032251                        {
    23042252                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
    2305                                 register int64 foo=(fNormal)[jj];
    2306                                 if( ( (*check)[jj]  == /*fNormal[jj]*/foo )
    2307                                  || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) )
    2308                                 {
    2309                                         notPar=TRUE;
    2310                                         break;
    2311                                 }
     2253//                              register int64 foo=(fNormal)[jj];
     2254//                              if( ( (*check)[jj]  == /*fNormal[jj]*/foo )
     2255//                               || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) )
     2256//                              {
     2257//                                      notPar=TRUE;
     2258//                                      break;
     2259//                              }
    23122260                        }
    23132261                        omFree(v);
    2314                         if (isParallel(*check,fNormal))Found a parallel vector. Add it
    2315                         if(notPar==FALSE)
    2316                         {
    2317                                 initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));pAdd = p_Add_q destroys args
     2262                        if (isParallel(*check,fNormal))//Found a parallel vector. Add it
     2263//                      if(notPar==FALSE)
     2264                        {
     2265                                initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));//pAdd = p_Add_q destroys args
    23182266                        }
    23192267                        delete check;
    2320                 }while
    2321                 omFree(v);
     2268                }//while
     2269//              omFree(v);
    23222270#ifdef gfan_DEBUG
    2323                 cout << "Initial Form=";
    2324                 pWrite(initialFormElement[ii]);
    2325                 cout << "---" << endl;
     2271//              cout << "Initial Form=";                               
     2272//              pWrite(initialFormElement[ii]);
     2273//              cout << "---" << endl;
    23262274#endif
    23272275                /*Now initialFormElement must be added to (ideal)initialForm */
     
    23292277                pDelete(&initialFormElement);
    23302278                omFree(leadExpV);
    2331         }for
     2279        }//for
    23322280#ifdef gfanp
    23332281        gettimeofday(&end, 0);
     
    23432291 * compute the factors \f$ a_i \f$
    23442292 */
    2345 NOTE: Should be replaced by kNF or kNF2
    2346 NOTE: Done
    2347 NOTE: removed with r12286
    2348 
     2293//NOTE: Should be replaced by kNF or kNF2
     2294//NOTE: Done
     2295//NOTE: removed with r12286
     2296               
    23492297/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
    23502298*/
    2351 NOTE: use kNF or kNF2 instead of restOfDivision
     2299//NOTE: use kNF or kNF2 instead of restOfDivision
    23522300inline ideal gcone::ffG(const ideal &H, const ideal &G)
    23532301{
     
    23562304        for (int ii=0;ii<size;ii++)
    23572305        {
    2358                 poly temp1;//=pInit();
    2359                 poly temp2;//=pInit();
    2360                 poly temp3;=pInit();//polys to temporarily store values for pSub
    2361                 res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
    2362                 temp1=pCopy(H->m[ii]);//TRY
    2363                 temp2=pCopy(res->m[ii]);
    2364                 NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
    2365                 temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));//TRY
    2366                 temp3=pSub(temp1, temp2);//TRY
    2367                 temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0)));NOTRY
     2306//              poly temp1;//=pInit();
     2307//              poly temp2;//=pInit();
     2308                poly temp3;//=pInit();//polys to temporarily store values for pSub
     2309//              res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
     2310//              temp1=pCopy(H->m[ii]);//TRY
     2311//              temp2=pCopy(res->m[ii]);
     2312                //NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
     2313//              temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));//TRY
     2314//              temp3=pSub(temp1, temp2);//TRY
     2315                temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0)));//NOTRY
    23682316                res->m[ii]=pCopy(temp3);
    2369                 res->m[ii]=pSub(temp1,temp2); //buggy
    2370                 cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
    2371                 pDelete(&temp1);//TRY
    2372                 pDelete(&temp2);
     2317                //res->m[ii]=pSub(temp1,temp2); //buggy         
     2318                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
     2319//              pDelete(&temp1);//TRY
     2320//              pDelete(&temp2);
    23732321                pDelete(&temp3);
    23742322        }
    23752323        return res;
    23762324}
    2377 
     2325       
    23782326/** \brief Preprocessing of inequlities
    23792327* Do some preprocessing on the matrix of inequalities
     
    23882336        int num_elts=0;
    23892337        int offset=0;*/
    2390         Remove zeroes (and strictly pos rows?)
     2338        //Remove zeroes (and strictly pos rows?)
    23912339        for(int ii=0;ii<ddineq->rowsize;ii++)
    23922340        {
    23932341                int64vec *iv = new int64vec(this->numVars);
    2394                 int64vec *ivNull = new int64vec(this->numVars);Needed for intvec64::compare(*int64vec)
     2342                int64vec *ivNull = new int64vec(this->numVars);//Needed for intvec64::compare(*int64vec)
    23952343                int posCtr=0;
    23962344                for(int jj=0;jj<this->numVars;jj++)
    23972345                {
    23982346                        (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
    2399                         if((*iv)[jj]>0)check for strictly pos rows
     2347                        if((*iv)[jj]>0)//check for strictly pos rows
    24002348                                posCtr++;
    2401                         Behold! This will delete the row for the standard simplex!
     2349                        //Behold! This will delete the row for the standard simplex!
    24022350                }
    2403                 if( (iv->compare(0)==0) || (posCtr==iv->length()) )
    2404                 if( (posCtr==iv->length()) || (iv->compare(ivNull)==0) )
     2351//              if( (iv->compare(0)==0) || (posCtr==iv->length()) )
     2352                if( (posCtr==iv->length()) || (iv->compare(ivNull)==0) )               
    24052353                {
    24062354                        dd_MatrixRowRemove(&ddineq,ii+1);
    2407                         ii--;Yes. This is on purpose
     2355                        ii--;//Yes. This is on purpose
    24082356                }
    24092357                delete iv;
    24102358                delete ivNull;
    24112359        }
    2412         Remove duplicates of rows
    2413         posRowsArray=NULL;
    2414         num_alloc=0;
    2415         num_elts=0;
    2416         offset=0;
    2417         int num_newRows = ddineq->rowsize;
    2418         for(int ii=0;ii<ddineq->rowsize-1;ii++)
    2419         for(int ii=0;ii<num_newRows-1;ii++)
    2420         {
    2421                 int64vec *iv = new int64vec(this->numVars);//1st vector to check against
    2422                 for(int jj=0;jj<this->numVars;jj++)
    2423                         (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
    2424                 for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++)
    2425                 {
    2426                         int64vec *ivCheck = new int64vec(this->numVars);//Checked against iv
    2427                         for(int kk=0;kk<this->numVars;kk++)
    2428                                 (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
    2429                         if (iv->compare(ivCheck)==0)
    2430                         {
    2431  //                             cout << "=" << endl;
    2432  //                             num_alloc++;
    2433  //                             void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int)));
    2434  //                             if(!tmp)
    2435  //                             {
    2436  //                                     WerrorS("Woah dude! Couldn't realloc memory\n");
    2437  //                                     exit(-1);
    2438  //                             }
    2439  //                             posRowsArray = (int*)tmp;
    2440  //                             posRowsArray[num_elts]=jj;
    2441  //                             num_elts++;
    2442                                 dd_MatrixRowRemove(&ddineq,jj+1);
    2443                                 num_newRows = ddineq->rowsize;
    2444                         }
    2445                         delete ivCheck;
    2446                 }
    2447                 delete iv;
    2448         }
    2449         for(int ii=0;ii<num_elts;ii++)
    2450         {
    2451                 dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
    2452                 offset++;
    2453         }
    2454         free(posRowsArray);
    2455         Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/
    2456 }preprocessInequalities
    2457 
     2360        //Remove duplicates of rows
     2361//      posRowsArray=NULL;
     2362//      num_alloc=0;
     2363//      num_elts=0;
     2364//      offset=0;
     2365//      int num_newRows = ddineq->rowsize;
     2366//      for(int ii=0;ii<ddineq->rowsize-1;ii++)
     2367//      for(int ii=0;ii<num_newRows-1;ii++)
     2368//      {
     2369//              int64vec *iv = new int64vec(this->numVars);//1st vector to check against
     2370//              for(int jj=0;jj<this->numVars;jj++)
     2371//                      (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
     2372//              for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++)
     2373//              {
     2374//                      int64vec *ivCheck = new int64vec(this->numVars);//Checked against iv
     2375//                      for(int kk=0;kk<this->numVars;kk++)
     2376//                              (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
     2377//                      if (iv->compare(ivCheck)==0)
     2378//                      {
     2379// //                           cout << "=" << endl;
     2380// //                           num_alloc++;
     2381// //                           void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int)));
     2382// //                           if(!tmp)
     2383// //                           {
     2384// //                                   WerrorS("Woah dude! Couldn't realloc memory\n");
     2385// //                                   exit(-1);
     2386// //                           }
     2387// //                           posRowsArray = (int*)tmp;
     2388// //                           posRowsArray[num_elts]=jj;
     2389// //                           num_elts++;
     2390//                              dd_MatrixRowRemove(&ddineq,jj+1);
     2391//                              num_newRows = ddineq->rowsize;
     2392//                      }
     2393//                      delete ivCheck;
     2394//              }
     2395//              delete iv;
     2396//      }
     2397//      for(int ii=0;ii<num_elts;ii++)
     2398//      {
     2399//              dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
     2400//              offset++;
     2401//      }
     2402//      free(posRowsArray);
     2403        //Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/       
     2404}//preprocessInequalities
     2405       
    24582406/** \brief Compute a Groebner Basis
    24592407 *
     
    24622410 *\return void
    24632411 */
    2464 inline void gcone::getGB(const ideal &inputIdeal)
     2412inline void gcone::getGB(const ideal &inputIdeal)               
    24652413{
    24662414        BITSET save=test;
     
    24712419        idNorm(gb);
    24722420        idSkipZeroes(gb);
    2473         this->gcBasis=gb;       write the GB into gcBasis
     2421        this->gcBasis=gb;       //write the GB into gcBasis
    24742422        test=save;
    2475 }void getGB
    2476 
     2423}//void getGB
     2424               
    24772425/** \brief Compute the negative of a given int64vec
    2478         */
     2426        */             
    24792427static int64vec* ivNeg(/*const*/int64vec *iv)
    2480 {       Hm, switching to int64vec const int64vec does no longer work
    2481         int64vec *res; = new int64vec(iv->length());
     2428{       //Hm, switching to int64vec const int64vec does no longer work
     2429        int64vec *res;// = new int64vec(iv->length());
    24822430        res=iv64Copy(iv);
    2483         *res *= (int)-1;
     2431        *res *= (int)-1;                                               
    24842432        return res;
    24852433}
     
    24892437*
    24902438*/
    2491 static int dotProduct(const int64vec &iva, const int64vec &ivb)
    2492 {
    2493         int res=0;
     2439static int dotProduct(const int64vec &iva, const int64vec &ivb)                         
     2440{                       
     2441        int res=0;     
    24942442        for (int i=0;i<pVariables;i++)
    24952443        {
    2496  #ifndef NDEBUG
    2497         (const_cast<int64vec*>(&iva))->show(1,0); (const_cast<int64vec*>(&ivb))->show(1,0);
    2498  #endif
     2444// #ifndef NDEBUG
     2445//      (const_cast<int64vec*>(&iva))->show(1,0); (const_cast<int64vec*>(&ivb))->show(1,0);
     2446// #endif
    24992447                res = res+(iva[i]*ivb[i]);
    25002448        }
     
    25062454 */
    25072455static bool isParallel(const int64vec &a,const int64vec &b)
    2508 {
    2509 /*#ifdef gfanp
     2456{       
     2457/*#ifdef gfanp 
    25102458        timeval start, end;
    25112459        gettimeofday(&start, 0);
    2512 #endif*/
     2460#endif*/               
    25132461        bool res;
    25142462        int lhs=dotProduct(a,b)*dotProduct(a,b);
    25152463        int rhs=dotProduct(a,a)*dotProduct(b,b);
    2516  #ifdef gfanp
    2517         gettimeofday(&end, 0);
    2518         t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    2519  #endif
    2520         return res;
     2464// #ifdef gfanp
     2465//      gettimeofday(&end, 0);
     2466//      t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
     2467// #endif       
     2468//      return res;
    25212469        return res = (lhs==rhs)?TRUE:FALSE;
    2522 }bool isParallel
     2470}//bool isParallel
    25232471
    25242472/** \brief Compute an interior point of a given cone
    2525  * Result will be written into int64vec iv.
     2473 * Result will be written into int64vec iv. 
    25262474 * Any rational point is automatically converted into an integer.
    25272475 */
    2528 void gcone::interiorPoint( dd_MatrixPtr &M, int64vec &iv) no const &M here since we want to remove redundant rows
     2476void gcone::interiorPoint( dd_MatrixPtr &M, int64vec &iv) //no const &M here since we want to remove redundant rows
    25292477{
    25302478        dd_LPPtr lp,lpInt;
     
    25322480        dd_LPSolverType solver=dd_DualSimplex;
    25332481        dd_LPSolutionPtr lpSol=NULL;
    2534         dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
    2535         dd_rowindex ddnewpos;
    2536         dd_NumberType numb;
    2537         M->representation=dd_Inequality;
    2538 
    2539         NOTE: Make this n-dimensional!
    2540         dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
    2541 
     2482//      dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
     2483//      dd_rowindex ddnewpos;
     2484        dd_NumberType numb;     
     2485        //M->representation=dd_Inequality;
     2486                       
     2487        //NOTE: Make this n-dimensional!
     2488        //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
     2489                                                                       
    25422490        /*NOTE: Leave the following line commented out!
    25432491        * Otherwise it will slow down computations a great deal
    25442492        * */
    2545         dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
    2546         if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
     2493//      dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
     2494        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
    25472495        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
    25482496        int jj=1;
     
    25502498        {
    25512499                dd_set_si(posRestr->matrix[ii][jj],1);
    2552                 jj++;
     2500                jj++;                                                   
    25532501        }
    25542502        dd_MatrixAppendTo(&M,posRestr);
     
    25582506        if (lp==NULL){WerrorS("LP is NULL");}
    25592507#ifdef gfan_DEBUG
    2560                         dd_WriteLP(stdout,lp);
    2561 #endif
    2562 
     2508//                      dd_WriteLP(stdout,lp);
     2509#endif 
     2510                                       
    25632511        lpInt=dd_MakeLPforInteriorFinding(lp);
    25642512        if (err!=dd_NoError){WerrorS("Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint");}
    25652513#ifdef gfan_DEBUG
    2566                         dd_WriteLP(stdout,lpInt);
    2567 #endif
    2568         dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
     2514//                      dd_WriteLP(stdout,lpInt);
     2515#endif                 
     2516//      dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
    25692517        if (err!=dd_NoError)
    25702518        {
    25712519                WerrorS("Error during dd_FindRelativeInterior in gcone::interiorPoint");
    25722520                dd_WriteErrorMessages(stdout, err);
    2573         }
    2574         dd_LPSolve(lpInt,solver,&err);  This will not result in a point from the relative interior
    2575         if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}
     2521        }                       
     2522        dd_LPSolve(lpInt,solver,&err);  //This will not result in a point from the relative interior
     2523//      if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}                                                                       
    25762524        lpSol=dd_CopyLPSolution(lpInt);
    2577         if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
     2525//      if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
    25782526#ifdef gfan_DEBUG
    25792527        printf("Interior point: ");
     
    25832531        }
    25842532        printf("\n");
    2585 #endif
    2586         NOTE The following strongly resembles parts of makeInt.
    2587         Maybe merge sometimes
     2533#endif                 
     2534        //NOTE The following strongly resembles parts of makeInt.
     2535        //Maybe merge sometimes
    25882536        mpz_t kgV; mpz_init(kgV);
    25892537        mpz_set_str(kgV,"1",10);
     
    26092557        }
    26102558#ifdef gfan_DEBUG
    2611                         iv.show();
    2612                         cout << endl;
     2559//                      iv.show();
     2560//                      cout << endl;
    26132561#endif
    26142562        mpq_clear(qkgV);
    26152563        mpz_clear(tmp);
    26162564        mpz_clear(den);
    2617         mpz_clear(kgV);
    2618 
     2565        mpz_clear(kgV);                 
     2566                       
    26192567        dd_FreeLPSolution(lpSol);
    26202568        dd_FreeLPData(lpInt);
    26212569        dd_FreeLPData(lp);
    2622         set_free(ddlinset);
    2623         set_free(ddredrows);
    2624 
    2625 }void interiorPoint(dd_MatrixPtr const &M)
     2570//      set_free(ddlinset);
     2571//      set_free(ddredrows);   
     2572                       
     2573}//void interiorPoint(dd_MatrixPtr const &M)
    26262574
    26272575/** Computes an interior point of a cone by taking two interior points a,b from two different facets
    26282576* and then computing b+(a-b)/2
    2629 * Of course this only works for flippable facets
    2630 * Two cases may occur:
     2577* Of course this only works for flippable facets 
     2578* Two cases may occur: 
    26312579* 1st: There are only two facets who share the only strictly positive ray
    26322580* 2nd: There are at least two facets which have a distinct positive ray
     
    26372585* positive => these lie in a plane and thus their sum is not from relative interior.
    26382586* So let's just sum up all rays, find one strictly positive and shift the point along that ray
    2639 *
     2587* 
    26402588* Used by noRevS
    26412589*NOTE no longer used nor maintained. MM Mar 9, 2010
    26422590*/
    2643  void gcone::interiorPoint2()
    2644  {//idPrint(this->gcBasis);
    2645  #ifdef gfan_DEBUG
    2646         if(this->ivIntPt!=NULL)
    2647                 WarnS("Interior point already exists - ovrewriting!");
    2648  #endif
    2649         facet *f1 = this->facetPtr;
    2650         facet *f2 = NULL;
    2651         int64vec *intF1=NULL;
    2652         while(f1!=NULL)
    2653         {
    2654                 if(f1->isFlippable)
    2655                 {
    2656                         facet *f1Ray = f1->codim2Ptr;
    2657                         while(f1Ray!=NULL)
    2658                         {
    2659                                 const int64vec *check=f1Ray->getRef2FacetNormal();
    2660                                 if(iv64isStrictlyPositive(check))
    2661                                 {
    2662                                         intF1=iv64Copy(check);
    2663                                         break;
    2664                                 }
    2665                                 f1Ray=f1Ray->next;
    2666                         }
    2667                 }
    2668                 if(intF1!=NULL)
    2669                         break;
    2670                 f1=f1->next;
    2671         }
    2672         if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
    2673                 f2=f1->next;
    2674         else
    2675                 f2=this->facetPtr;
    2676         if(intF1==NULL && hasHomInput==TRUE)
    2677         {
    2678                 intF1 = new int64vec(this->numVars);
    2679                 for(int ii=0;ii<this->numVars;ii++)
    2680                         (*intF1)[ii]=1;
    2681         }
    2682         assert(f1); assert(f2);
    2683         int64vec *intF2=f2->getInteriorPoint();
    2684         mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
    2685         mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
    2686         mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
    2687         for(int ii=0;ii<this->numVars;ii++)
    2688         {
    2689                 mpq_init(qPosRay[ii]);
    2690                 mpq_init(qIntPt[ii]);
    2691                 mpq_init(qPosIntPt[ii]);
    2692         }
    2693         //Compute a+((b-a)/2) && Convert intF1 to mpq
    2694         for(int ii=0;ii<this->numVars;ii++)
    2695         {
    2696                 mpq_t a,b;
    2697                 mpq_init(a); mpq_init(b);
    2698                 mpq_set_si(a,(*intF1)[ii],1);
    2699                 mpq_set_si(b,(*intF2)[ii],1);
    2700                 mpq_t diff;
    2701                 mpq_init(diff);
    2702                 mpq_sub(diff,b,a);      //diff=b-a
    2703                 mpq_t quot;
    2704                 mpq_init(quot);
    2705                 mpq_div_2exp(quot,diff,1);      //quot=diff/2=(b-a)/2
    2706                 mpq_clear(diff);
    2707                 //Don't be clever and reuse diff here
    2708                 mpq_t sum; mpq_init(sum);
    2709                 mpq_add(sum,b,quot);    //sum=b+quot=a+(b-a)/2
    2710                 mpq_set(qIntPt[ii],sum);
    2711                 mpq_clear(sum);
    2712                 mpq_clear(quot);
    2713                 mpq_clear(a); mpq_clear(b);
    2714                 //Now for intF1
    2715                 mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
    2716         }
    2717         //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
    2718         while(TRUE)
    2719         {
    2720                 bool success=FALSE;
    2721                 int posCtr=0;
    2722                 for(int ii=0;ii<this->numVars;ii++)
    2723                 {
    2724                         mpq_t sum; mpq_init(sum);
    2725                         mpq_add(sum,qPosRay[ii],qIntPt[ii]);
    2726                         mpq_set(qPosIntPt[ii],sum);
    2727                         mpq_clear(sum);
    2728                         if(mpq_sgn(qPosIntPt[ii])==1)
    2729                                 posCtr++;
    2730                 }
    2731                 if(posCtr==this->numVars)//qPosIntPt > 0
    2732                         break;
    2733                 else
    2734                 {
    2735                         mpq_t qTwo; mpq_init(qTwo);
    2736                         mpq_set_ui(qTwo,2,1);
    2737                         for(int jj=0;jj<this->numVars;jj++)
    2738                         {
    2739                                 mpq_t tmp; mpq_init(tmp);
    2740                                 mpq_mul(tmp,qPosRay[jj],qTwo);
    2741                                 mpq_set( qPosRay[jj], tmp);
    2742                                 mpq_clear(tmp);
    2743                         }
    2744                         mpq_clear(qTwo);
    2745                 }
    2746         }//while
    2747         //Now qPosIntPt ought to be >0, so convert back to int :D
    2748         /*Compute lcm of the denominators*/
    2749         mpz_t *denom = new mpz_t[this->numVars];
    2750         mpz_t tmp,kgV;
    2751         mpz_init(tmp); mpz_init(kgV);
    2752         for (int ii=0;ii<this->numVars;ii++)
    2753         {
    2754                 mpz_t z;
    2755                 mpz_init(z);
    2756                 mpq_get_den(z,qPosIntPt[ii]);
    2757                 mpz_init(denom[ii]);
    2758                 mpz_set( denom[ii], z);
    2759                 mpz_clear(z);
    2760         }
    2761 
    2762         mpz_set(tmp,denom[0]);
    2763         for (int ii=0;ii<this->numVars;ii++)
    2764         {
    2765                 mpz_lcm(kgV,tmp,denom[ii]);
    2766                 mpz_set(tmp,kgV);
    2767         }
    2768         mpz_clear(tmp);
    2769         /*Multiply the nominators by kgV*/
    2770         mpq_t qkgV,res;
    2771         mpq_init(qkgV);
    2772         mpq_canonicalize(qkgV);
    2773         mpq_init(res);
    2774         mpq_canonicalize(res);
    2775 
    2776         mpq_set_num(qkgV,kgV);
    2777         int64vec *n=new int64vec(this->numVars);
    2778         for (int ii=0;ii<this->numVars;ii++)
    2779         {
    2780                 mpq_canonicalize(qPosIntPt[ii]);
    2781                 mpq_mul(res,qkgV,qPosIntPt[ii]);
    2782                 (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
    2783         }
    2784         this->setIntPoint(n);
    2785         delete n;
    2786         delete [] qPosIntPt;
    2787         delete [] denom;
    2788         delete [] qPosRay;
    2789         delete [] qIntPt;
    2790         mpz_clear(kgV);
    2791         mpq_clear(qkgV); mpq_clear(res);
    2792  }
    2793 
     2591// void gcone::interiorPoint2()
     2592// {//idPrint(this->gcBasis);
     2593// #ifdef gfan_DEBUG
     2594//      if(this->ivIntPt!=NULL)
     2595//              WarnS("Interior point already exists - ovrewriting!");
     2596// #endif
     2597//      facet *f1 = this->facetPtr;
     2598//      facet *f2 = NULL;
     2599//      int64vec *intF1=NULL;
     2600//      while(f1!=NULL)
     2601//      {
     2602//              if(f1->isFlippable)
     2603//              {
     2604//                      facet *f1Ray = f1->codim2Ptr;
     2605//                      while(f1Ray!=NULL)
     2606//                      {
     2607//                              const int64vec *check=f1Ray->getRef2FacetNormal();
     2608//                              if(iv64isStrictlyPositive(check))
     2609//                              {
     2610//                                      intF1=iv64Copy(check);
     2611//                                      break;
     2612//                              }                               
     2613//                              f1Ray=f1Ray->next;
     2614//                      }
     2615//              }
     2616//              if(intF1!=NULL)
     2617//                      break;
     2618//              f1=f1->next;
     2619//      }
     2620//      if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
     2621//              f2=f1->next;
     2622//      else
     2623//              f2=this->facetPtr;
     2624//      if(intF1==NULL && hasHomInput==TRUE)
     2625//      {
     2626//              intF1 = new int64vec(this->numVars);
     2627//              for(int ii=0;ii<this->numVars;ii++)
     2628//                      (*intF1)[ii]=1;
     2629//      }
     2630//      assert(f1); assert(f2);
     2631//      int64vec *intF2=f2->getInteriorPoint();
     2632//      mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
     2633//      mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
     2634//      mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
     2635//      for(int ii=0;ii<this->numVars;ii++)
     2636//      {
     2637//              mpq_init(qPosRay[ii]);
     2638//              mpq_init(qIntPt[ii]);
     2639//              mpq_init(qPosIntPt[ii]);
     2640//      }       
     2641//      //Compute a+((b-a)/2) && Convert intF1 to mpq
     2642//      for(int ii=0;ii<this->numVars;ii++)
     2643//      {
     2644//              mpq_t a,b;
     2645//              mpq_init(a); mpq_init(b);
     2646//              mpq_set_si(a,(*intF1)[ii],1);
     2647//              mpq_set_si(b,(*intF2)[ii],1);
     2648//              mpq_t diff;
     2649//              mpq_init(diff);
     2650//              mpq_sub(diff,b,a);      //diff=b-a
     2651//              mpq_t quot;
     2652//              mpq_init(quot);
     2653//              mpq_div_2exp(quot,diff,1);      //quot=diff/2=(b-a)/2
     2654//              mpq_clear(diff);
     2655//              //Don't be clever and reuse diff here
     2656//              mpq_t sum; mpq_init(sum);
     2657//              mpq_add(sum,b,quot);    //sum=b+quot=a+(b-a)/2
     2658//              mpq_set(qIntPt[ii],sum);
     2659//              mpq_clear(sum);
     2660//              mpq_clear(quot);
     2661//              mpq_clear(a); mpq_clear(b);
     2662//              //Now for intF1
     2663//              mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
     2664//      }
     2665//      //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
     2666//      while(TRUE)
     2667//      {       
     2668//              bool success=FALSE;
     2669//              int posCtr=0;   
     2670//              for(int ii=0;ii<this->numVars;ii++)
     2671//              {
     2672//                      mpq_t sum; mpq_init(sum);
     2673//                      mpq_add(sum,qPosRay[ii],qIntPt[ii]);
     2674//                      mpq_set(qPosIntPt[ii],sum);
     2675//                      mpq_clear(sum);
     2676//                      if(mpq_sgn(qPosIntPt[ii])==1)
     2677//                              posCtr++;
     2678//              }
     2679//              if(posCtr==this->numVars)//qPosIntPt > 0
     2680//                      break;
     2681//              else
     2682//              {
     2683//                      mpq_t qTwo; mpq_init(qTwo);
     2684//                      mpq_set_ui(qTwo,2,1);
     2685//                      for(int jj=0;jj<this->numVars;jj++)
     2686//                      {
     2687//                              mpq_t tmp; mpq_init(tmp);
     2688//                              mpq_mul(tmp,qPosRay[jj],qTwo);
     2689//                              mpq_set( qPosRay[jj], tmp);
     2690//                              mpq_clear(tmp);
     2691//                      }
     2692//                      mpq_clear(qTwo);
     2693//              }
     2694//      }//while
     2695//      //Now qPosIntPt ought to be >0, so convert back to int :D
     2696//      /*Compute lcm of the denominators*/
     2697//      mpz_t *denom = new mpz_t[this->numVars];
     2698//      mpz_t tmp,kgV;
     2699//      mpz_init(tmp); mpz_init(kgV);
     2700//      for (int ii=0;ii<this->numVars;ii++)
     2701//      {
     2702//              mpz_t z;
     2703//              mpz_init(z);
     2704//              mpq_get_den(z,qPosIntPt[ii]);
     2705//              mpz_init(denom[ii]);
     2706//              mpz_set( denom[ii], z);
     2707//              mpz_clear(z);                           
     2708//      }
     2709//             
     2710//      mpz_set(tmp,denom[0]);
     2711//      for (int ii=0;ii<this->numVars;ii++)
     2712//      {
     2713//              mpz_lcm(kgV,tmp,denom[ii]);
     2714//              mpz_set(tmp,kgV);                               
     2715//      }
     2716//      mpz_clear(tmp);
     2717//      /*Multiply the nominators by kgV*/
     2718//      mpq_t qkgV,res;
     2719//      mpq_init(qkgV);
     2720//      mpq_canonicalize(qkgV);         
     2721//      mpq_init(res);
     2722//      mpq_canonicalize(res);
     2723//                             
     2724//      mpq_set_num(qkgV,kgV);
     2725//      int64vec *n=new int64vec(this->numVars);
     2726//      for (int ii=0;ii<this->numVars;ii++)
     2727//      {
     2728//              mpq_canonicalize(qPosIntPt[ii]);
     2729//              mpq_mul(res,qkgV,qPosIntPt[ii]);
     2730//              (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
     2731//      }
     2732//      this->setIntPoint(n);
     2733//      delete n;
     2734//      delete [] qPosIntPt;
     2735//      delete [] denom;
     2736//      delete [] qPosRay;
     2737//      delete [] qIntPt;
     2738//      mpz_clear(kgV);
     2739//      mpq_clear(qkgV); mpq_clear(res);
     2740// }
     2741       
    27942742/** \brief Copy a ring and add a weighted ordering in first place
    2795  *
     2743 * 
    27962744 */
    2797 ring gcone::rCopyAndAddWeight(const ring &r, int64vec *ivw)
     2745ring gcone::rCopyAndAddWeight(const ring &r, int64vec *ivw)                             
    27982746{
    27992747        ring res=rCopy0(r);
    28002748        int jj;
    2801 
     2749                       
    28022750        omFree(res->order);
    28032751        res->order =(int *)omAlloc0(4*sizeof(int/*64*/));
     
    28082756        omfree(res->wvhdl);
    28092757        res->wvhdl =(int **)omAlloc0(4*sizeof(int/*64*/**));
    2810 
     2758                       
    28112759        res->order[0]=ringorder_a/*64*/;
    28122760        res->block0[0]=1;
    28132761        res->block1[0]=res->N;
    2814         res->order[1]=ringorder_dp;     basically useless, since that should never be used
     2762        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                   
    28152763        res->block0[1]=1;
    28162764        res->block1[1]=res->N;
     
    28202768        int/*64*/ *A=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
    28212769        for (jj=0;jj<length;jj++)
    2822         {
     2770        {                               
    28232771                A[jj]=(*ivw)[jj];
    28242772                if((*ivw)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::rCopyAndAddWeight!\n");
    2825         }
     2773        }                       
    28262774        res->wvhdl[0]=(int *)A;
    28272775        res->block1[0]=length;
    2828 
     2776                       
    28292777        rComplete(res);
    28302778        return res;
    2831 }rCopyAndAdd
    2832 
    2833 ring gcone::rCopyAndAddWeight2(const ring &r,const int64vec *ivw, const int64vec *fNormal)
     2779}//rCopyAndAdd
     2780               
     2781ring gcone::rCopyAndAddWeight2(const ring &r,const int64vec *ivw, const int64vec *fNormal)                             
    28342782{
    28352783        ring res=rCopy0(r);
    2836 
     2784                       
    28372785        omFree(res->order);
    28382786        res->order =(int *)omAlloc0(5*sizeof(int/*64*/));
     
    28432791        omfree(res->wvhdl);
    28442792        res->wvhdl =(int **)omAlloc0(5*sizeof(int/*64*/**));
    2845 
     2793                       
    28462794        res->order[0]=ringorder_a/*64*/;
    28472795        res->block0[0]=1;
     
    28502798        res->block0[1]=1;
    28512799        res->block1[1]=res->N;
    2852 
     2800       
    28532801        res->order[2]=ringorder_dp;
    28542802        res->block0[2]=1;
    28552803        res->block1[2]=res->N;
    2856 
     2804       
    28572805        res->order[3]=ringorder_C;
    28582806
     
    28612809        int/*64*/ *A2=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
    28622810        for (int jj=0;jj<length;jj++)
    2863         {
     2811        {                               
    28642812                A1[jj]=(*ivw)[jj];
    28652813                A2[jj]=-(*fNormal)[jj];
    28662814                if((*ivw)[jj]>=INT_MAX || (*fNormal)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::rCopyAndAddWeight2!\n");
    2867         }
     2815        }                       
    28682816        res->wvhdl[0]=(int *)A1;
    28692817        res->block1[0]=length;
     
    28732821        return res;
    28742822}
    2875 
    2876 NOTE not needed anywhere
    2877  ring rCopyAndChangeWeight(ring const &r, int64vec *ivw)
    2878  {
    2879         ring res=rCopy0(currRing);
    2880         rComplete(res);
    2881         rSetWeightVec(res,(int64*)ivw);
    2882         //rChangeCurrRing(rnew);
    2883         return res;
    2884  }
    2885 
     2823               
     2824//NOTE not needed anywhere
     2825// ring rCopyAndChangeWeight(ring const &r, int64vec *ivw)
     2826// {
     2827//      ring res=rCopy0(currRing);
     2828//      rComplete(res);
     2829//      rSetWeightVec(res,(int64*)ivw);
     2830//      //rChangeCurrRing(rnew);
     2831//      return res;
     2832// }
     2833               
    28862834/** \brief Checks whether a given facet is a search facet
    28872835 * Determines whether a given facet of a cone is the search facet of a neighbouring cone
     
    28902838 * that is first crossed during the generic walk.
    28912839 * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
    2892  * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
     2840 * If this is the case, then our facet is indeed a search facet and TRUE is retuned. 
    28932841*/
    2894 removed with r12286
    2895 
     2842//removed with r12286
     2843               
    28962844/** \brief Check for equality of two intvecs
    28972845 */
     
    29092857        return res;
    29102858}
    2911 
     2859               
    29122860/** \brief The reverse search algorithm
    29132861 */
    2914 removed with r12286
     2862//removed with r12286
    29152863/** \brief Compute the lineality/homogeneity space
    29162864* It is the kernel of the inequality matrix Ax=0
     
    29212869        dd_MatrixPtr res;
    29222870        dd_MatrixPtr ddineq;
    2923         ddineq=dd_CopyMatrix(this->ddFacets);
    2924         Add a row of 0s in 0th place
     2871        ddineq=dd_CopyMatrix(this->ddFacets);   
     2872        //Add a row of 0s in 0th place
    29252873        dd_MatrixPtr ddAppendRowOfZeroes=dd_CreateMatrix(1,this->numVars+1);
    29262874        dd_MatrixPtr ddFoo=dd_AppendMatrix(ddAppendRowOfZeroes,ddineq);
     
    29292877        ddineq=dd_CopyMatrix(ddFoo);
    29302878        dd_FreeMatrix(ddFoo);
    2931         Cohen starts here
    2932         int dimKer=0;Cohen calls this r
    2933         int m=ddineq->rowsize;Rows
    2934         int n=ddineq->colsize;Cols
     2879        //Cohen starts here
     2880        int dimKer=0;//Cohen calls this r
     2881        int m=ddineq->rowsize;//Rows
     2882        int n=ddineq->colsize;//Cols
    29352883        int c[m+1];
    29362884        int d[n+1];
     
    29382886                c[ii]=0;
    29392887        for(int ii=0;ii<n;ii++)
    2940                 d[ii]=0;
    2941 
     2888                d[ii]=0;       
     2889       
    29422890        for(int k=1;k<n;k++)
    29432891        {
    2944                 //Let's find a j s.t. m[j][k]!=0 && c[j]=0
    2945                 int condCtr=0;Check each row for zeroness
     2892                //Let's find a j s.t. m[j][k]!=0 && c[j]=0             
     2893                int condCtr=0;//Check each row for zeroness
    29462894                for(int j=1;j<m;j++)
    29472895                {
     
    29552903                                mpq_set(ratd,quot);
    29562904                                mpq_canonicalize(ratd);
    2957 
     2905               
    29582906                                mpq_set_str(ddineq->matrix[j][k],"-1",10);
    29592907                                for(int ss=k+1;ss<n;ss++)
    29602908                                {
    29612909                                        mpq_t prod; mpq_init(prod);
    2962                                         mpq_mul(prod, ratd, ddineq->matrix[j][ss]);
     2910                                        mpq_mul(prod, ratd, ddineq->matrix[j][ss]);                             
    29632911                                        mpq_set(ddineq->matrix[j][ss],prod);
    29642912                                        mpq_canonicalize(ddineq->matrix[j][ss]);
    29652913                                        mpq_clear(prod);
    2966                                 }
     2914                                }               
    29672915                                for(int ii=1;ii<m;ii++)
    29682916                                {
     
    29702918                                        {
    29712919                                                mpq_set(ratd,ddineq->matrix[ii][k]);
    2972                                                 mpq_set_str(ddineq->matrix[ii][k],"0",10);
     2920                                                mpq_set_str(ddineq->matrix[ii][k],"0",10);                     
    29732921                                                for(int ss=k+1;ss<n;ss++)
    29742922                                                {
     
    29842932                                        }
    29852933                                }
    2986                                 c[j]=k;
     2934                                c[j]=k;         
    29872935                                d[k]=j;
    2988                                 mpq_clear(quot); mpq_clear(ratd); mpq_clear(one);
     2936                                mpq_clear(quot); mpq_clear(ratd); mpq_clear(one);       
    29892937                        }
    29902938                        else
    29912939                                condCtr++;
    2992                 }
    2993                 if(condCtr==m-1)        Why -1 ???
     2940                }                       
     2941                if(condCtr==m-1)        //Why -1 ???
    29942942                {
    29952943                        dimKer++;
    29962944                        d[k]=0;
    2997                         break;//goto _4;
    2998                 }else{
     2945//                      break;//goto _4;
     2946                }//else{
    29992947                /*Eliminate*/
    3000         }for(k
    3001         /*Output kernel, i.e. set gcone::dd_LinealitySpace here*/
    3002         gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
     2948        }//for(k
     2949        /*Output kernel, i.e. set gcone::dd_LinealitySpace here*/       
     2950//      gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
    30032951        res = dd_CreateMatrix(dimKer,this->numVars+1);
    30042952        int row=-1;
     
    30122960                                if(d[ii]>0)
    30132961                                        mpq_set(res->matrix[row][ii],ddineq->matrix[d[ii]][kk]);
    3014                                 else if(ii==kk)
     2962                                else if(ii==kk)                         
    30152963                                        mpq_set_str(res->matrix[row][ii],"1",10);
    30162964                                mpq_canonicalize(res->matrix[row][ii]);
     
    30202968        dd_FreeMatrix(ddineq);
    30212969        return(res);
    3022         Better safe than sorry:
    3023         NOTE dd_LinealitySpace contains RATIONAL ENTRIES
    3024         Therefore if you need integer ones: CALL gcone::makeInt(...) method
    3025 }
     2970        //Better safe than sorry:
     2971        //NOTE dd_LinealitySpace contains RATIONAL ENTRIES
     2972        //Therefore if you need integer ones: CALL gcone::makeInt(...) method
     2973}       
    30262974
    30272975
     
    30342982 */
    30352983void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
    3036 {
    3037         facet *SearchListRoot = new facet(); The list containing ALL facets we come across
     2984{       
     2985        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
    30382986        facet *SearchListAct;
    30392987        SearchListAct = NULL;
    3040         SearchListAct = SearchListRoot;
     2988        //SearchListAct = SearchListRoot;                       
    30412989        gcone *gcAct;
    30422990        gcAct = &gcRoot;
    3043         gcone *gcPtr;   Pointer to end of linked list of cones
     2991        gcone *gcPtr;   //Pointer to end of linked list of cones
    30442992        gcPtr = &gcRoot;
    3045         gcone *gcNext;  Pointer to next cone to be visited
     2993        gcone *gcNext;  //Pointer to next cone to be visited
    30462994        gcNext = &gcRoot;
    30472995        gcone *gcHead;
    3048         gcHead = &gcRoot;
     2996        gcHead = &gcRoot;                       
    30492997        facet *fAct;
    3050         fAct = gcAct->facetPtr;
     2998        fAct = gcAct->facetPtr;                         
    30512999        ring rAct;
    30523000        rAct = currRing;
    3053         int UCNcounter=gcAct->getUCN();
     3001        int UCNcounter=gcAct->getUCN(); 
    30543002#ifdef gfan_DEBUG
    30553003        printf("NoRevs\n");
    30563004        printf("Facets are:\n");
    30573005        gcAct->showFacets();
    3058 #endif
     3006#endif                 
    30593007        /*Compute lineality space here
    30603008        * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/
    30613009        if(dd_LinealitySpace==NULL)
    30623010                dd_LinealitySpace = gcAct->computeLinealitySpace();
    3063         /*End of lineality space computation*/
    3064         gcAct->getCodim2Normals(*gcAct);
     3011        /*End of lineality space computation*/         
     3012        //gcAct->getCodim2Normals(*gcAct);
    30653013        if(fAct->codim2Ptr==NULL)
    3066                 gcAct->getExtremalRays(*gcAct);
     3014                gcAct->getExtremalRays(*gcAct);         
    30673015        /* Make a copy of the facet list of first cone
    30683016           Since the operations getCodim2Normals and normalize affect the facets
    3069            we must not memcpy them before these ops! */
    3070         /*facet *codim2Act; codim2Act = NULL;
     3017           we must not memcpy them before these ops! */ 
     3018        /*facet *codim2Act; codim2Act = NULL;                   
    30713019        facet *sl2Root;
    3072         facet *sl2Act;*/
     3020        facet *sl2Act;*/                       
    30733021        for(int ii=0;ii<this->numFacets;ii++)
    30743022        {
    3075                 only copy flippable facets! NOTE: We do not need flipRing or any such information.
     3023                //only copy flippable facets! NOTE: We do not need flipRing or any such information.
    30763024                if(fAct->isFlippable==TRUE)
    30773025                {
    30783026                        /*Using shallow copy here*/
    30793027#ifdef SHALLOW
    3080                         if( ii==0 || (ii>0 && SearchListAct==NULL) ) 1st facet may be non-flippable
     3028                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
    30813029                        {
    30823030                                if(SearchListRoot!=NULL) delete(SearchListRoot);
    30833031                                SearchListRoot = fAct->shallowCopy(*fAct);
    3084                                 SearchListAct = SearchListRoot; SearchListRoot is already 'new'ed at beginning of method!
     3032                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
    30853033                        }
    30863034                        else
    3087                         {
     3035                        {                       
    30883036                                SearchListAct->next = fAct->shallowCopy(*fAct);
    3089                                 SearchListAct = SearchListAct->next;
     3037                                SearchListAct = SearchListAct->next;                                           
    30903038                        }
    30913039                        fAct=fAct->next;
     
    30943042                        int64vec *fNormal;
    30953043                        fNormal = fAct->getFacetNormal();
    3096                         if( ii==0 || (ii>0 && SearchListAct==NULL) ) 1st facet may be non-flippable
    3097                         {
    3098                                 SearchListAct = SearchListRoot; SearchListRoot is already 'new'ed at beginning of method!
     3044                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
     3045                        {                                               
     3046                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
    30993047                        }
    31003048                        else
    3101                         {
     3049                        {                       
    31023050                                SearchListAct->next = new facet();
    3103                                 SearchListAct = SearchListAct->next;
     3051                                SearchListAct = SearchListAct->next;                                           
    31043052                        }
    31053053                        SearchListAct->setFacetNormal(fNormal);
     
    31073055                        SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
    31083056                        SearchListAct->isFlippable=TRUE;
    3109                         Copy int point as well
     3057                        //Copy int point as well
    31103058                        int64vec *ivIntPt;
    31113059                        ivIntPt = fAct->getInteriorPoint();
    31123060                        SearchListAct->setInteriorPoint(ivIntPt);
    31133061                        delete(ivIntPt);
    3114                         Copy codim2-facets
     3062                        //Copy codim2-facets
    31153063                        facet *codim2Act;
    3116                         codim2Act = NULL;
     3064                        codim2Act = NULL;                       
    31173065                        facet *sl2Root;
    3118                         facet *sl2Act;
     3066                        facet *sl2Act;                 
    31193067                        codim2Act=fAct->codim2Ptr;
    31203068                        SearchListAct->codim2Ptr = new facet(2);
    31213069                        sl2Root = SearchListAct->codim2Ptr;
    3122                         sl2Act = sl2Root;
     3070                        sl2Act = sl2Root;                       
    31233071                        for(int jj=0;jj<fAct->numCodim2Facets;jj++)
    3124                         for(int jj=0;jj<fAct->numRays-1;jj++)
     3072//                      for(int jj=0;jj<fAct->numRays-1;jj++)
    31253073                        {
    31263074                                int64vec *f2Normal;
    31273075                                f2Normal = codim2Act->getFacetNormal();
    31283076                                if(jj==0)
    3129                                 {
     3077                                {                                               
    31303078                                        sl2Act = sl2Root;
    31313079                                        sl2Act->setFacetNormal(f2Normal);
     
    31373085                                        sl2Act->setFacetNormal(f2Normal);
    31383086                                }
    3139                                 delete f2Normal;
     3087                                delete f2Normal;                               
    31403088                                codim2Act = codim2Act->next;
    31413089                        }
     
    31433091                        delete fNormal;
    31443092#endif
    3145                 }if(fAct->isFlippable==TRUE)
     3093                }//if(fAct->isFlippable==TRUE)                 
    31463094                else {fAct = fAct->next;}
    3147         }End of copying facets into SLA
    3148 
    3149         SearchListAct = SearchListRoot; Set to beginning of list
     3095        }//End of copying facets into SLA
     3096       
     3097        SearchListAct = SearchListRoot; //Set to beginning of list
    31503098        /*Make SearchList doubly linked*/
    31513099#ifndef NDEBUG
     
    31643112                if(SearchListAct->next!=NULL){
    31653113#endif
    3166                         SearchListAct->next->prev = SearchListAct;
     3114                        SearchListAct->next->prev = SearchListAct;                                     
    31673115                }
    31683116                SearchListAct = SearchListAct->next;
    31693117        }
    3170         SearchListAct = SearchListRoot; Set to beginning of List
    3171 
    3172         fAct = gcAct->facetPtr;//???
    3173         gcAct->writeConeToFile(*gcAct);
     3118        SearchListAct = SearchListRoot; //Set to beginning of List
     3119       
     3120//      fAct = gcAct->facetPtr;//???
     3121        gcAct->writeConeToFile(*gcAct);                 
    31743122        /*End of initialisation*/
    3175 
     3123       
    31763124        fAct = SearchListAct;
    31773125        /*2nd step
    31783126          Choose a facet from SearchList, flip it and forget the previous cone
    31793127          We always choose the first facet from SearchList as facet to be flipped
    3180         */
    3181         while( (SearchListAct!=NULL))&& counter<490)
    3182         {NOTE See to it that the cone is only changed after ALL facets have been flipped!
    3183                 fAct = SearchListAct;
     3128        */     
     3129        while( (SearchListAct!=NULL))//&& counter<490)
     3130        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
     3131                fAct = SearchListAct;           
    31843132                while(fAct!=NULL)
    3185                 while( (fAct->getUCN() == fAct->next->getUCN()) )
    3186                 {       Since SLA should only contain flippables there should be no need to check for that
     3133//              while( (fAct->getUCN() == fAct->next->getUCN()) )               
     3134                {       //Since SLA should only contain flippables there should be no need to check for that                   
    31873135                        gcAct->flip2(gcAct->gcBasis,fAct);
    3188                         NOTE rCopy needed?
     3136                        //NOTE rCopy needed?
    31893137                        ring rTmp=rCopy(fAct->flipRing);
    31903138                        rComplete(rTmp);
    31913139                        rChangeCurrRing(rTmp);
    3192                         gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);copy constructor!
     3140                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);//copy constructor!
    31933141                        /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
    31943142                         * Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
     
    32003148                        */
    32013149                        idDelete((ideal *)&fAct->flipGB);
    3202                         rDelete(fAct->flipRing);
     3150                        rDelete(fAct->flipRing);                       
    32033151                        gcTmp->getConeNormals(gcTmp->gcBasis/*, FALSE*/);
    3204                         gcTmp->getCodim2Normals(*gcTmp);
     3152//                      gcTmp->getCodim2Normals(*gcTmp);
    32053153                        gcTmp->getExtremalRays(*gcTmp);
    3206                         NOTE Order rays lex here
    3207                         gcTmp->orderRays();
    3208                         //NOTE If flip2 is used we need to get an interior point of gcTmp
    3209                         // and replace gcTmp->baseRing with an appropriate ring with only
    3210                         // one weight
    3211                         gcTmp->interiorPoint2();
     3154                        //NOTE Order rays lex here
     3155                        gcTmp->orderRays();                     
     3156//                      //NOTE If flip2 is used we need to get an interior point of gcTmp
     3157//                      // and replace gcTmp->baseRing with an appropriate ring with only
     3158//                      // one weight
     3159//                      gcTmp->interiorPoint2();
    32123160                        gcTmp->replaceDouble_ringorder_a_ByASingleOne();
    3213                         gcTmp->ddFacets should not be needed anymore, so
     3161                        //gcTmp->ddFacets should not be needed anymore, so
    32143162                        dd_FreeMatrix(gcTmp->ddFacets);
    32153163#ifdef gfan_DEBUG
    3216                         gcTmp->showFacets(1);
     3164//                      gcTmp->showFacets(1);
    32173165#endif
    32183166                        /*add facets to SLA here*/
     
    32283176  #endif
    32293177                        SearchListRoot=tmp;
    3230                         SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
    3231 #else
     3178                        //SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
     3179#else 
    32323180                        SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
    3233 #endif ifdef SHALLOW
    3234                         gcTmp->writeConeToFile(*gcTmp);
     3181#endif //ifdef SHALLOW
     3182//                      gcTmp->writeConeToFile(*gcTmp);
    32353183                        if(gfanHeuristic==1)
    32363184                        {
    32373185                                gcTmp->writeConeToFile(*gcTmp);
    3238                                 idDelete((ideal*)&gcTmp->gcBasis);Whonder why?
     3186                                idDelete((ideal*)&gcTmp->gcBasis);//Whonder why?
    32393187                                rDelete(gcTmp->baseRing);
    3240                         }
     3188                        }                       
    32413189#ifdef gfan_DEBUG
    32423190                        if(SearchListRoot!=NULL)
    32433191                                showSLA(*SearchListRoot);
    3244 #endif
     3192#endif                 
    32453193                        rChangeCurrRing(gcAct->baseRing);
    32463194                        rDelete(rTmp);
    3247                         doubly linked for easier removal
     3195                        //doubly linked for easier removal
    32483196                        gcTmp->prev = gcPtr;
    32493197                        gcPtr->next=gcTmp;
    32503198                        gcPtr=gcPtr->next;
    3251                         Cleverly disguised exit condition follows
     3199                        //Cleverly disguised exit condition follows
    32523200                        if(fAct->getUCN() == fAct->next->getUCN())
    32533201                        {
    32543202                                printf("Switching UCN from %i to %i\n",fAct->getUCN(),fAct->next->getUCN());
    3255                                 fAct=fAct->next;
     3203                                fAct=fAct->next;                               
    32563204                        }
    32573205                        else
    32583206                        {
    3259                                 rDelete(gcAct->baseRing);
    3260                                 printf("break\n");
     3207                                //rDelete(gcAct->baseRing);
     3208//                              printf("break\n");
    32613209                                break;
    32623210                        }
    3263                         fAct=fAct->next;
    3264                 }while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
    3265                 Search for cone with smallest UCN
     3211//                      fAct=fAct->next;
     3212                }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );         
     3213                //Search for cone with smallest UCN
    32663214#ifndef NDEBUG
    3267   #if SIZEOF_LONG==8    64 bit
     3215  #if SIZEOF_LONG==8    //64 bit
    32683216                while(gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && SearchListRoot!=NULL)
    32693217  #elif SIZEOF_LONG == 4
     
    32723220#endif
    32733221#ifdef NDEBUG
    3274                 while(gcNext!=NULL && SearchListRoot!=NULL)
    3275 #endif
    3276                 {
     3222                while(gcNext!=NULL && SearchListRoot!=NULL)     
     3223#endif
     3224                {                       
    32773225                        if( gcNext->getUCN() == SearchListRoot->getUCN() )
    32783226                        {
     
    32803228                                {
    32813229                                        gcAct = gcNext;
    3282                                         Seems better to not to use rCopy here
    3283                                         rAct=rCopy(gcAct->baseRing);
     3230                                        //Seems better to not to use rCopy here
     3231//                                      rAct=rCopy(gcAct->baseRing);
    32843232                                        rAct=gcAct->baseRing;
    32853233                                        rComplete(rAct);
    3286                                         rChangeCurrRing(rAct);
     3234                                        rChangeCurrRing(rAct);                                         
    32873235                                        break;
    32883236                                }
     
    32903238                                {
    32913239                                        gcone *gcDel;
    3292                                         gcDel = gcAct;
     3240                                        gcDel = gcAct;                                 
    32933241                                        gcAct = gcNext;
    3294                                         Read st00f from file &
    3295                                         implant the GB into gcAct
     3242                                        //Read st00f from file &
     3243                                        //implant the GB into gcAct
    32963244                                        readConeFromFile(gcAct->getUCN(), gcAct);
    3297                                         Kill the baseRing but ONLY if it is not the ring the computation started in!
    3298                                         if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
    3299                                                 rDelete(gcDel->baseRing);
    3300                                         rAct=rCopy(gcAct->baseRing);
    3301                                         /*The ring change occurs in readConeFromFile, so as to
     3245                                        //Kill the baseRing but ONLY if it is not the ring the computation started in!
     3246//                                      if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
     3247//                                              rDelete(gcDel->baseRing);
     3248//                                      rAct=rCopy(gcAct->baseRing);
     3249                                        /*The ring change occurs in readConeFromFile, so as to 
    33023250                                        assure that gcAct->gcBasis belongs to the right ring*/
    33033251                                        rAct=gcAct->baseRing;
     
    33053253                                        rChangeCurrRing(rAct);
    33063254                                        break;
    3307                                 }
     3255                                }                               
    33083256                        }
    33093257                        else if(gcNext->getUCN() < SearchListRoot->getUCN() )
    33103258                        {
    3311                                 idDelete( (ideal*)&gcNext->gcBasis );
    3312                                 rDelete(gcNext->baseRing);//TODO Why does this crash?
     3259                                idDelete( (ideal*)&gcNext->gcBasis );                           
     3260//                              rDelete(gcNext->baseRing);//TODO Why does this crash?
    33133261                        }
    33143262                        /*else
     
    33203268                                        if(gcDel->getUCN()!=1)
    33213269                                                rDelete(gcDel->baseRing);
    3322                                 }
    3323                         }*/
     3270                                }                                       
     3271                        }*/                     
    33243272                        gcNext = gcNext->next;
    33253273                }
    33263274                UCNcounter++;
    3327                 SearchListAct = SearchListRoot;
     3275                SearchListAct = SearchListRoot;         
    33283276        }
    33293277        printf("\nFound %i cones - terminating\n", counter);
    3330 }void noRevS(gcone &gc)
    3331 
    3332 
     3278}//void noRevS(gcone &gc)       
     3279               
     3280               
    33333281/** \brief Make a set of rational vectors into integers
    33343282 *
     
    33393287 */
    33403288void gcone::makeInt(const dd_MatrixPtr &M, const int line, int64vec &n)
    3341 {
     3289{                       
    33423290        mpz_t *denom = new mpz_t[this->numVars];
    33433291        for(int ii=0;ii<this->numVars;ii++)
     
    33453293                mpz_init_set_str(denom[ii],"0",10);
    33463294        }
    3347 
     3295               
    33483296        mpz_t kgV,tmp;
    33493297        mpz_init(kgV);
    33503298        mpz_init(tmp);
    3351 
     3299                       
    33523300        for (int ii=0;ii<(M->colsize)-1;ii++)
    33533301        {
     
    33563304                mpq_get_den(z,M->matrix[line-1][ii+1]);
    33573305                mpz_set( denom[ii], z);
    3358                 mpz_clear(z);
    3359         }
    3360 
     3306                mpz_clear(z);                           
     3307        }
     3308                                               
    33613309        /*Compute lcm of the denominators*/
    33623310        mpz_set(tmp,denom[0]);
     
    33643312        {
    33653313                mpz_lcm(kgV,tmp,denom[ii]);
    3366                 mpz_set(tmp,kgV);
    3367         }
    3368         mpz_clear(tmp);
     3314                mpz_set(tmp,kgV);                               
     3315        }
     3316        mpz_clear(tmp); 
    33693317        /*Multiply the nominators by kgV*/
    33703318        mpq_t qkgV,res;
     
    33723320        mpq_set_str(qkgV,"1",10);
    33733321        mpq_canonicalize(qkgV);
    3374 
     3322                       
    33753323        mpq_init(res);
    33763324        mpq_set_str(res,"1",10);
    33773325        mpq_canonicalize(res);
    3378 
     3326                       
    33793327        mpq_set_num(qkgV,kgV);
    3380 
    3381         mpq_canonicalize(qkgV);
    3382         int ggT=1;
     3328                       
     3329//      mpq_canonicalize(qkgV);
     3330//      int ggT=1;
    33833331        for (int ii=0;ii<(M->colsize)-1;ii++)
    33843332        {
    33853333                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
    33863334                n[ii]=(int64)mpz_get_d(mpq_numref(res));
    3387                 ggT=int64gcd(ggT,n[ii]);
     3335//              ggT=int64gcd(ggT,n[ii]);
    33883336        }
    33893337        int64 ggT=n[0];
    33903338        for(int ii=0;ii<this->numVars;ii++)
    3391                 ggT=int64gcd(ggT,n[ii]);
    3392         Normalisation
     3339                ggT=int64gcd(ggT,n[ii]);       
     3340        //Normalisation
    33933341        if(ggT>1)
    33943342        {
     
    33983346        delete [] denom;
    33993347        mpz_clear(kgV);
    3400         mpq_clear(qkgV); mpq_clear(res);
    3401 
     3348        mpq_clear(qkgV); mpq_clear(res);                       
     3349                       
    34023350}
    34033351/**
    3404  * For all codim-2-facets we compute the gcd of the components of the facet normal and
     3352 * For all codim-2-facets we compute the gcd of the components of the facet normal and 
    34053353 * divide it out. Thus we get a normalized representation of each
    34063354 * (codim-2)-facet normal, i.e. a primitive vector
    34073355 * Actually we now also normalize the facet normals.
    34083356 */
    3409  void gcone::normalize()
    3410  {
    3411         int *ggT = new int;
    3412                 *ggT=1;
    3413         facet *fAct;
    3414         facet *codim2Act;
    3415         fAct = this->facetPtr;
    3416         codim2Act = fAct->codim2Ptr;
    3417         while(fAct!=NULL)
    3418         {
    3419                 int64vec *fNormal;
    3420                 fNormal = fAct->getFacetNormal();
    3421                 int *ggT = new int;
    3422                 *ggT=1;
    3423                 for(int ii=0;ii<this->numVars;ii++)
    3424                 {
    3425                         *ggT=intgcd((*ggT),(*fNormal)[ii]);
    3426                 }
    3427                 if(*ggT>1)//We only need to do this if the ggT is non-trivial
    3428                 {
    3429  //                     int64vec *fCopy = fAct->getFacetNormal();
    3430                         for(int ii=0;ii<this->numVars;ii++)
    3431                                 (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT);
    3432                         fAct->setFacetNormal(fNormal);
    3433                 }
    3434                 delete fNormal;
    3435                 delete ggT;
    3436                 /*And now for the codim2*/
    3437                 while(codim2Act!=NULL)
    3438                 {
    3439                         int64vec *n;
    3440                         n=codim2Act->getFacetNormal();
    3441                         int *ggT=new int;
    3442                         *ggT=1;
    3443                         for(int ii=0;ii<this->numVars;ii++)
    3444                         {
    3445                                 *ggT = intgcd((*ggT),(*n)[ii]);
    3446                         }
    3447                         if(*ggT>1)
    3448                         {
    3449                                 for(int ii=0;ii<this->numVars;ii++)
    3450                                 {
    3451                                         (*n)[ii] = ((*n)[ii])/(*ggT);
    3452                                 }
    3453                                 codim2Act->setFacetNormal(n);
    3454                         }
    3455                         codim2Act = codim2Act->next;
    3456                         delete n;
    3457                         delete ggT;
    3458                 }
    3459                 fAct = fAct->next;
    3460         }
    3461  }
    3462 
    3463 /** \brief Enqueue new facets into the searchlist
     3357// void gcone::normalize()
     3358// {
     3359//      int *ggT = new int;
     3360//              *ggT=1;
     3361//      facet *fAct;
     3362//      facet *codim2Act;
     3363//      fAct = this->facetPtr;
     3364//      codim2Act = fAct->codim2Ptr;
     3365//      while(fAct!=NULL)
     3366//      {
     3367//              int64vec *fNormal;
     3368//              fNormal = fAct->getFacetNormal();
     3369//              int *ggT = new int;
     3370//              *ggT=1;
     3371//              for(int ii=0;ii<this->numVars;ii++)
     3372//              {
     3373//                      *ggT=intgcd((*ggT),(*fNormal)[ii]);
     3374//              }
     3375//              if(*ggT>1)//We only need to do this if the ggT is non-trivial
     3376//              {
     3377// //                   int64vec *fCopy = fAct->getFacetNormal();
     3378//                      for(int ii=0;ii<this->numVars;ii++)
     3379//                              (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT);
     3380//                      fAct->setFacetNormal(fNormal);
     3381//              }               
     3382//              delete fNormal;
     3383//              delete ggT;
     3384//              /*And now for the codim2*/
     3385//              while(codim2Act!=NULL)
     3386//              {                               
     3387//                      int64vec *n;
     3388//                      n=codim2Act->getFacetNormal();
     3389//                      int *ggT=new int;
     3390//                      *ggT=1;
     3391//                      for(int ii=0;ii<this->numVars;ii++)
     3392//                      {
     3393//                              *ggT = intgcd((*ggT),(*n)[ii]);
     3394//                      }
     3395//                      if(*ggT>1)
     3396//                      {
     3397//                              for(int ii=0;ii<this->numVars;ii++)
     3398//                              {
     3399//                                      (*n)[ii] = ((*n)[ii])/(*ggT);
     3400//                              }
     3401//                              codim2Act->setFacetNormal(n);
     3402//                      }
     3403//                      codim2Act = codim2Act->next;
     3404//                      delete n;
     3405//                      delete ggT;
     3406//              }
     3407//              fAct = fAct->next;
     3408//      }
     3409// }
     3410
     3411/** \brief Enqueue new facets into the searchlist 
    34643412 * The searchlist (SLA for short) is implemented as a FIFO
    34653413 * Checks are done as follows:
    34663414 * 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
    3467  * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
     3415 * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the 
    34683416 *      facet from SLA and do not add fAct.
    34693417 * It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
     
    34803428        facet *slHead;
    34813429        slHead = f;
    3482         facet *slAct;   called with f=SearchListRoot
     3430        facet *slAct;   //called with f=SearchListRoot
    34833431        slAct = f;
    3484         facet *slEnd;   Pointer to end of SLA
     3432        facet *slEnd;   //Pointer to end of SLA
    34853433        slEnd = f;
    3486         facet *slEndStatic;     //marks the end before a new facet is added
     3434//      facet *slEndStatic;     //marks the end before a new facet is added             
    34873435        facet *fAct;
    34883436        fAct = this->facetPtr;
     
    34943442        volatile facet *deleteMarker;
    34953443        deleteMarker = NULL;
    3496 
     3444                       
    34973445        /** \brief  Flag to mark a facet that might be added
    34983446         * The following case may occur:
     
    35003448         * This does however not mean that there does not exist a facet behind the current slAct that is actually equal
    35013449         * to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
    3502          * FALSE and the according slAct is deleted.
     3450         * FALSE and the according slAct is deleted. 
    35033451         * If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
    35043452         */
    35053453
    3506         /**A facet was removed, lengthOfSearchlist-- thus we must not rely on
     3454        /**A facet was removed, lengthOfSearchlist-- thus we must not rely on 
    35073455         * if(notParallelCtr==lengthOfSearchList) but rather
    35083456         * if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
     
    35133461                slEnd=slEnd->next;
    35143462        }
    3515         /*1st step: compare facetNormals*/
     3463        /*1st step: compare facetNormals*/                     
    35163464        while(fAct!=NULL)
    3517         {
     3465        {                                               
    35183466                if(fAct->isFlippable==TRUE)
    35193467                {
    35203468                        int64vec *fNormal=NULL;
    3521                         fNormal=fAct->getFacetNormal();
     3469                        fNormal=fAct->getFacetNormal();                 
    35223470                        slAct = slHead;
    3523                         /*If slAct==NULL and fAct!=NULL
     3471                        /*If slAct==NULL and fAct!=NULL 
    35243472                        we just copy all remaining facets into SLA*/
    35253473                        if(slAct==NULL)
     
    35293477                                while(fCopy!=NULL)
    35303478                                {
    3531                                         if(fCopy->isFlippable==TRUE)We must assure fCopy is also flippable
     3479                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
    35323480                                        {
    35333481                                                if(slAct==NULL)
    35343482                                                {
    3535                                                         slAct = new facet(*fCopy/*,TRUE*/);copy constructor
    3536                                                         slHead = slAct;
     3483                                                        slAct = new facet(*fCopy/*,TRUE*/);//copy constructor
     3484                                                        slHead = slAct;                                                         
    35373485                                                }
    35383486                                                else
     
    35443492                                        fCopy = fCopy->next;
    35453493                                }
    3546                                 break;Where does this lead to?
     3494                                break;//Where does this lead to?
    35473495                        }
    35483496                        /*End of dumping into SLA*/
     
    35543502#ifdef gfan_DEBUG
    35553503                                printf("Checking facet (");fNormal->show(1,1);printf(") against (");slNormal->show(1,1);printf(")\n");
    3556 #endif
    3557                                 if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
    3558                                         exit(-1);
     3504#endif                         
     3505//                              if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
     3506//                                      exit(-1);
    35593507                                if(areEqual2(fAct,slAct))
    3560                                 {
     3508                                {                                       
    35613509                                        deleteMarker = slAct;
    35623510                                        if(slAct==slHead)
    3563                                         {
    3564                                                 slHead = slAct->next;
     3511                                        {                                               
     3512                                                slHead = slAct->next;                                           
    35653513                                                if(slHead!=NULL)
    35663514                                                        slHead->prev = NULL;
     
    35703518                                                slEnd=slEnd->prev;
    35713519                                                slEnd->next = NULL;
    3572                                         }
     3520                                        }                                                               
    35733521                                        else
    35743522                                        {
     
    35803528                                        if(deleteMarker!=NULL)
    35813529                                        {
    3582                                                 delete deleteMarker;
    3583                                                 deleteMarker=NULL;
     3530//                                              delete deleteMarker;
     3531//                                              deleteMarker=NULL;
    35843532                                        }
    35853533#ifdef gfan_DEBUG
     
    35873535#endif
    35883536                                        delete slNormal;
    3589                                         break;leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
     3537                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
    35903538                                }
    35913539                                slAct = slAct->next;
     
    35963544                                slAct = slAct->next;*/
    35973545                                if(deleteMarker!=NULL)
    3598                                 {
    3599                                         delete deleteMarker;
    3600                                         deleteMarker=NULL;
     3546                                {                                               
     3547//                                      delete deleteMarker;
     3548//                                      deleteMarker=NULL;
    36013549                                }
    36023550                                delete slNormal;
    3603                                                 if slAct was marked as to be deleted, delete it here!
    3604                         }while(slAct!=NULL)
     3551                                                //if slAct was marked as to be deleted, delete it here!
     3552                        }//while(slAct!=NULL)
    36053553                        if(removalOccured==FALSE)
    36063554                        {
    36073555#ifdef gfan_DEBUG
    3608                                 cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl;
    3609 #endif
    3610                                 Add fAct to SLA
     3556//                              cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl;
     3557#endif
     3558                                //Add fAct to SLA
    36113559                                facet *marker;
    36123560                                marker = slEnd;
     
    36153563
    36163564                                slEnd->next = new facet();
    3617                                 slEnd = slEnd->next;Update slEnd
     3565                                slEnd = slEnd->next;//Update slEnd
    36183566                                facet *slEndCodim2Root;
    36193567                                facet *slEndCodim2Act;
    36203568                                slEndCodim2Root = slEnd->codim2Ptr;
    36213569                                slEndCodim2Act = slEndCodim2Root;
    3622 
     3570                                                               
    36233571                                slEnd->setUCN(this->getUCN());
    36243572                                slEnd->isFlippable = TRUE;
    36253573                                slEnd->setFacetNormal(fNormal);
    3626                                 NOTE Add interior point here
    3627                                 This is ugly but needed for flip2
    3628                                 Better: have slEnd->interiorPoint point to fAct->interiorPoint
    3629                                 NOTE Only reference -> c.f. copy constructor
    3630                                 slEnd->setInteriorPoint(fAct->getInteriorPoint());
     3574                                //NOTE Add interior point here
     3575                                //This is ugly but needed for flip2
     3576                                //Better: have slEnd->interiorPoint point to fAct->interiorPoint
     3577                                //NOTE Only reference -> c.f. copy constructor
     3578                                //slEnd->setInteriorPoint(fAct->getInteriorPoint());
    36313579                                slEnd->interiorPoint = fAct->interiorPoint;
    36323580                                slEnd->prev = marker;
    3633                                 Copy codim2-facets
    3634                                 int64vec *f2Normal=new int64vec(this->numVars);
     3581                                //Copy codim2-facets
     3582                                //int64vec *f2Normal=new int64vec(this->numVars);
    36353583                                while(f2Act!=NULL)
    36363584                                {
     
    36403588                                        {
    36413589                                                slEndCodim2Root = new facet(2);
    3642                                                 slEnd->codim2Ptr = slEndCodim2Root;
     3590                                                slEnd->codim2Ptr = slEndCodim2Root;                     
    36433591                                                slEndCodim2Root->setFacetNormal(f2Normal);
    36443592                                                slEndCodim2Act = slEndCodim2Root;
    36453593                                        }
    3646                                         else
     3594                                        else                                   
    36473595                                        {
    36483596                                                slEndCodim2Act->next = new facet(2);
     
    36533601                                        delete f2Normal;
    36543602                                }
    3655                                 gcone::lengthOfSearchList++;
    3656                         }if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
     3603                                gcone::lengthOfSearchList++;                                                   
     3604                        }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
    36573605                        fAct = fAct->next;
    36583606                        delete fNormal;
    3659                         delete slNormal;
    3660                 }if(fAct->isFlippable==TRUE)
     3607//                      delete slNormal;
     3608                }//if(fAct->isFlippable==TRUE)
    36613609                else
    36623610                {
     
    36653613                if(gcone::maxSize<gcone::lengthOfSearchList)
    36663614                        gcone::maxSize= gcone::lengthOfSearchList;
    3667         }while(fAct!=NULL)
     3615        }//while(fAct!=NULL)
    36683616#ifdef gfanp
    36693617        gettimeofday(&end, 0);
    36703618        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    3671 #endif
     3619#endif                                         
    36723620        return slHead;
    3673 }addC2N
     3621}//addC2N
    36743622
    36753623/** Enqueuing using shallow copies*/
     
    36833631        facet *slHead;
    36843632        slHead = f;
    3685         facet *slAct;   called with f=SearchListRoot
     3633        facet *slAct;   //called with f=SearchListRoot
    36863634        slAct = f;
    3687         static facet *slEnd;    Pointer to end of SLA
     3635        static facet *slEnd;    //Pointer to end of SLA
    36883636        if(slEnd==NULL)
    36893637                slEnd = f;
    3690 
     3638       
    36913639        facet *fAct;
    3692         fAct = this->facetPtr;New facets to compare
     3640        fAct = this->facetPtr;//New facets to compare
    36933641        facet *codim2Act;
    36943642        codim2Act = this->facetPtr->codim2Ptr;
     
    36973645        {
    36983646                slEnd=slEnd->next;
    3699         }
     3647        }       
    37003648        while(fAct!=NULL)
    37013649        {
     
    37083656                                printf("Zero length SLA\n");
    37093657                                facet *fCopy;
    3710                                 fCopy = fAct;
     3658                                fCopy = fAct;                           
    37113659                                while(fCopy!=NULL)
    37123660                                {
    3713                                         if(fCopy->isFlippable==TRUE)We must assure fCopy is also flippable
     3661                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
    37143662                                        {
    37153663                                                if(slAct==NULL)
    37163664                                                {
    3717                                                         slAct = slAct->shallowCopy(*fCopy);shallow copy constructor
    3718                                                         slHead = slAct;
     3665                                                        slAct = slAct->shallowCopy(*fCopy);//shallow copy constructor
     3666                                                        slHead = slAct;                                                         
    37193667                                                }
    37203668                                                else
     
    37263674                                        fCopy = fCopy->next;
    37273675                                }
    3728                                 break;  WTF?
     3676                                break;  //WTF?
    37293677                        }
    37303678                        /*Comparison starts here*/
     
    37343682#ifdef gfan_DEBUG
    37353683        printf("Checking facet (");fAct->fNormal->show(1,1);printf(") against (");slAct->fNormal->show(1,1);printf(")\n");
    3736 #endif
     3684#endif 
    37373685                                if(areEqual2(fAct,slAct))
    37383686                                {
     
    37403688                                        if(slAct==slHead)
    37413689                                        {
    3742                                                 fDeleteMarker=slHead;
    3743                                                 printf("headUCN@enq=%i\n",slHead->getUCN());
    3744                                                 slHead = slAct->next;
    3745                                                 printf("headUCN@enq=%i\n",slHead->getUCN());
     3690//                                              fDeleteMarker=slHead;
     3691//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
     3692                                                slHead = slAct->next;                                           
     3693//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
    37463694                                                if(slHead!=NULL)
    37473695                                                {
     
    37493697                                                }
    37503698                                                fDeleteMarker->shallowDelete();
    3751                                                 delete fDeleteMarker;//NOTE this messes up fAct in noRevS!
    3752                                                 printf("headUCN@enq=%i\n",slHead->getUCN());
     3699                                                //delete fDeleteMarker;//NOTE this messes up fAct in noRevS!
     3700//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
    37533701                                        }
    37543702                                        else if (slAct==slEnd)
     
    37583706                                                fDeleteMarker->shallowDelete();
    37593707                                                delete(fDeleteMarker);
    3760                                         }
     3708                                        }                                                               
    37613709                                        else
    37623710                                        {
     
    37713719printf("Removing (");fAct->fNormal->show(1,1);printf(") from list\n");
    37723720#endif
    3773                                         break;leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
     3721                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
    37743722                                }
    3775                                 slAct = slAct->next;
    3776                         }while(slAct!=NULL)
     3723                                slAct = slAct->next;                           
     3724                        }//while(slAct!=NULL)
    37773725                        if(removalOccured==FALSE)
    37783726                        {
     
    37843732                        }
    37853733                        fAct = fAct->next;
    3786                         if(fDeleteMarker!=NULL)
    3787                         {
    3788                                 fDeleteMarker->shallowDelete();
    3789                                 delete(fDeleteMarker);
    3790                                 fDeleteMarker=NULL;
    3791                         }
     3734//                      if(fDeleteMarker!=NULL)
     3735//                      {
     3736//                              fDeleteMarker->shallowDelete();
     3737//                              delete(fDeleteMarker);
     3738//                              fDeleteMarker=NULL;
     3739//                      }
    37923740                }
    37933741                else
    37943742                        fAct = fAct->next;
    3795         }while(fAct!=NULL)
    3796 
     3743        }//while(fAct!=NULL)
     3744       
    37973745#ifdef gfanp
    37983746        gettimeofday(&end, 0);
    37993747        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    3800 #endif
    3801         printf("headUCN@enq=%i\n",slHead->getUCN());
     3748#endif 
     3749//      printf("headUCN@enq=%i\n",slHead->getUCN());
    38023750        return slHead;
    38033751}
    38043752
    3805 /**
     3753/** 
    38063754* During flip2 every gc->baseRing gets two ringorder_a
    38073755* To avoid having an awful lot of those in the end we endow
     
    38243772        omfree(replacementRing->wvhdl);
    38253773        replacementRing->wvhdl =(int **)omAlloc0(4*sizeof(int/*64*/**));
    3826 
     3774                       
    38273775        replacementRing->order[0]=ringorder_a/*64*/;
    38283776        replacementRing->block0[0]=1;
    38293777        replacementRing->block1[0]=replacementRing->N;
    3830 
     3778               
    38313779        replacementRing->order[1]=ringorder_dp;
    38323780        replacementRing->block0[1]=1;
    38333781        replacementRing->block1[1]=replacementRing->N;
    3834 
     3782       
    38353783        replacementRing->order[2]=ringorder_C;
    38363784
    3837         int64vec *ivw = this->getIntPoint(TRUE);returns a reference
    3838         assert(this->ivIntPt);
    3839         int length=ivw->length();
     3785        int64vec *ivw = this->getIntPoint(TRUE);//returns a reference
     3786//      assert(this->ivIntPt); 
     3787        int length=ivw->length();       
    38403788        int/*64*/ *A=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
    38413789        for (int jj=0;jj<length;jj++)
     
    38433791                A[jj]=(*ivw)[jj];
    38443792                if((*ivw)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::replaceDouble_ringorder_a_ByASingleOne!\n");
    3845         }
    3846         delete ivw; //Not needed if this->getIntPoint(TRUE)
     3793        }       
     3794        //delete ivw; //Not needed if this->getIntPoint(TRUE)
    38473795        replacementRing->wvhdl[0]=(int *)A;
    38483796        replacementRing->block1[0]=length;
     
    38533801        rDelete(this->baseRing);
    38543802        this->baseRing=rCopy(replacementRing);
    3855         this->gcBasis=idCopy(temporaryGroebnerBasis);
     3803        this->gcBasis=idCopy(temporaryGroebnerBasis);   
    38563804        /*And back to where we came from*/
    38573805        rChangeCurrRing(srcRing);
    38583806        idDelete( (ideal*)&temporaryGroebnerBasis );
    3859         rDelete(replacementRing);
     3807        rDelete(replacementRing);       
    38603808}
    38613809
     
    38773825        return p;
    38783826}
    3879 
     3827       
     3828static int intgcd(const int &a, const int &b)
     3829{
     3830        int r, p=a, q=b;
     3831        if(p < 0)
     3832                p = -p;
     3833        if(q < 0)
     3834                q = -q;
     3835        while(q != 0)
     3836        {
     3837                r = p % q;
     3838                p = q;
     3839                q = r;
     3840        }
     3841        return p;
     3842}
     3843               
    38803844/** \brief Construct a dd_MatrixPtr from a cone's list of facets
    38813845 * NO LONGER USED
    38823846 */
    3883  inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
    3884  {
    3885         facet *fAct;
    3886         fAct = gc.facetPtr;
    3887 
    3888         dd_MatrixPtr M;
    3889         dd_rowrange ddrows;
    3890         dd_colrange ddcols;
    3891         ddcols=(this->numVars)+1;
    3892         ddrows=this->numFacets;
    3893         dd_NumberType numb = dd_Integer;
    3894         M=dd_CreateMatrix(ddrows,ddcols);
    3895 
    3896         int jj=0;
    3897 
    3898         while(fAct!=NULL)
    3899         {
    3900                 int64vec *comp;
    3901                 comp = fAct->getFacetNormal();
    3902                 for(int ii=0;ii<this->numVars;ii++)
    3903                 {
    3904                         dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
    3905                 }
    3906                 jj++;
    3907                 delete comp;
    3908                 fAct=fAct->next;
    3909         }
    3910         return M;
    3911  }
    3912 
     3847// inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
     3848// {
     3849//      facet *fAct;
     3850//      fAct = gc.facetPtr;
     3851//     
     3852//      dd_MatrixPtr M;
     3853//      dd_rowrange ddrows;
     3854//      dd_colrange ddcols;
     3855//      ddcols=(this->numVars)+1;
     3856//      ddrows=this->numFacets;
     3857//      dd_NumberType numb = dd_Integer;
     3858//      M=dd_CreateMatrix(ddrows,ddcols);                       
     3859//                     
     3860//      int jj=0;
     3861//     
     3862//      while(fAct!=NULL)
     3863//      {
     3864//              int64vec *comp;
     3865//              comp = fAct->getFacetNormal();
     3866//              for(int ii=0;ii<this->numVars;ii++)
     3867//              {
     3868//                      dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
     3869//              }
     3870//              jj++;
     3871//              delete comp;
     3872//              fAct=fAct->next;                               
     3873//      }                       
     3874//      return M;
     3875// }
     3876               
    39133877/** \brief Write information about a cone into a file on disk
    39143878 *
     
    39263890        stringstream ss;
    39273891        ss << UCN;
    3928         string UCNstr = ss.str();
    3929 
     3892        string UCNstr = ss.str();               
     3893                       
    39303894        string prefix="/tmp/Singular/cone";
    3931         string prefix="./";     //crude hack -> run tests in separate directories
     3895//      string prefix="./";     //crude hack -> run tests in separate directories
    39323896        string suffix=".sg";
    39333897        string filename;
     
    39353899        filename.append(UCNstr);
    39363900        filename.append(suffix);
    3937 
    3938         int thisPID = getpid();
    3939         ss << thisPID;
    3940         string strPID = ss.str();
    3941         string prefix="./";
    3942 
     3901               
     3902//      int thisPID = getpid();
     3903//      ss << thisPID;
     3904//      string strPID = ss.str();
     3905//      string prefix="./";
     3906                                               
    39433907        ofstream gcOutputFile(filename.c_str());
    39443908        assert(gcOutputFile);
    39453909        facet *fAct;
    3946         fAct = gc.facetPtr;
     3910        fAct = gc.facetPtr;                     
    39473911        facet *f2Act;
    39483912        f2Act=fAct->codim2Ptr;
    3949 
     3913       
    39503914        char *ringString = rString(gc.baseRing);
    3951 
     3915                       
    39523916        if (!gcOutputFile)
    39533917        {
     
    39553919        }
    39563920        else
    3957         {
     3921        {       
    39583922                gcOutputFile << "UCN" << endl;
    39593923                gcOutputFile << this->UCN << endl;
    3960                 gcOutputFile << "RING" << endl;
     3924                gcOutputFile << "RING" << endl; 
    39613925                gcOutputFile << ringString << endl;
    39623926                gcOutputFile << "GCBASISLENGTH" << endl;
    3963                 gcOutputFile << IDELEMS(gc.gcBasis) << endl;
    3964                 Write this->gcBasis into file
    3965                 gcOutputFile << "GCBASIS" << endl;
     3927                gcOutputFile << IDELEMS(gc.gcBasis) << endl;                   
     3928                //Write this->gcBasis into file
     3929                gcOutputFile << "GCBASIS" << endl;                             
    39663930                for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
    3967                 {
     3931                {                                       
    39683932                        gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
    3969                 }
    3970 
    3971                 gcOutputFile << "FACETS" << endl;
    3972 
     3933                }                               
     3934                                       
     3935                gcOutputFile << "FACETS" << endl;                                                               
     3936               
    39733937                while(fAct!=NULL)
    3974                 {
     3938                {       
    39753939                        const int64vec *iv=fAct->getRef2FacetNormal();
    3976                         iv=fAct->getRef2FacetNormal();//->getFacetNormal();