source: git/kernel/gfan.cc @ 1f637e

spielwiese
Last change on this file since 1f637e was 1f637e, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: replaced pVariables by currRing->N... ADD: added 'void rChangeCurrRing(ring r);' to polys/polys.h CHG: started changing gr_kstd2.cc and syz*.cc
  • Property mode set to 100644
File size: 125.4 KB
RevLine 
[b3e45f]1/*
[6623f3]2Compute the Groebner fan of an ideal
[832b70]3$Author: monerjan $
[7c0ec6]4$Date: 2009/11/03 06:57:32 $
5$Header: /usr/local/Singular/cvsroot/kernel/gfan.cc,v 1.103 2009/11/03 06:57:32 monerjan Exp $
[341696]6$Id$
[b3e45f]7*/
8
[599326]9#include <kernel/mod2.h>
[34fcf8]10
[4199b3]11#ifdef HAVE_FANS
[0f401f]12#include <misc/options.h>
[599326]13#include <kernel/kstd1.h>
14#include <kernel/kutil.h>
[210e07]15#include <polys/polys.h>
[599326]16#include <kernel/ideals.h>
17#include <kernel/kmatrix.h>
18#include <kernel/GMPrat.h>
[4199b3]19
[76cfef]20#include <polys/monomials/ring.h>       //apparently not needed
21#include lists.h
22#include <polys/prCopy.h>
[599326]23#include <kernel/stairc.h>
[4199b3]24#include <fstream>      //read-write cones to files
[ca02ab]25#include <string>
26#include <sstream>
[5ece85]27#include <stdlib.h>
[be7ab3f]28#include <assert.h>
[4199b3]29#include <gfanlib/gfanlib.h>
[f85223]30
[1fcefd]31/*DO NOT REMOVE THIS*/
[f85223]32#ifndef GMPRATIONAL
33#define GMPRATIONAL
34#endif
[fc4f9a]35
[4199b3]36//Hacks for different working places
[8d82b4]37#define p800
[fc4f9a]38
[fc1971]39#ifdef p800
[b9a447]40#include <../cddlib-094f/lib-src-gmp/setoper.h>
41#include <../cddlib-094f/lib-src-gmp/cdd.h>
42#include <../cddlib-094f/lib-src-gmp/cddmp.h>
[fc1971]43#endif
44
[b3e45f]45#ifndef gfan_DEBUG
[4199b3]46// #define gfan_DEBUG
[fc1971]47#ifndef gfan_DEBUGLEVEL
48#define gfan_DEBUGLEVEL 1
49#endif
[b3e45f]50#endif
[2f72c55]51
[4199b3]52//NOTE Defining this will slow things down!
53//Only good for very coarse profiling
54// #define gfanp
[2f72c55]55#ifdef gfanp
[96e408]56  #include <sys/time.h>
57  #include <iostream>
[2f72c55]58#endif
59
[96e408]60//NOTE DO NOT REMOVE THIS
[73809b]61#ifndef SHALLOW
[96e408]62  #define SHALLOW
[73809b]63#endif
64
[4199b3]65#ifndef USE_ZFAN
[96e408]66  #define USE_ZFAN
[4199b3]67#endif
68
[ac6e45]69#include <gfan.h>
[1fcefd]70using namespace std;
[9f9b142]71
[147167f]72#define ivIsStrictlyPositive iv64isStrictlyPositive
73
[fe919c]74/**
75*\brief Class facet
[53e33d9]76*       Implements the facet structure as a linked list
[fe919c]77*
78*/
[4199b3]79               
[bb503c7]80/** \brief The default constructor for facets
81*/
[4199b3]82facet::facet()                 
[ac6e45]83{
[4199b3]84        // Pointer to next facet.  */
[ac6e45]85        /* Defaults to NULL. This way there is no need to check explicitly */
86        this->fNormal=NULL;
[4199b3]87        this->interiorPoint=NULL;               
[ac6e45]88        this->UCN=0;
89        this->codim2Ptr=NULL;
[4199b3]90        this->codim=1;          //default for (codim-1)-facets
[ac6e45]91        this->numCodim2Facets=0;
[5ece85]92        this->numRays=0;
[ac6e45]93        this->flipGB=NULL;
94        this->next=NULL;
[4199b3]95        this->prev=NULL;               
96        this->flipRing=NULL;    //the ring on the other side
[ac6e45]97        this->isFlippable=FALSE;
98}
[4199b3]99               
[ec6c52]100/** \brief Constructor for facets of codim == 2
[005d00a]101* Note that as of now the code of the constructors is only for facets and codim2-faces. One
102* could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
[ac6e45]103*/
[73809b]104facet::facet(const int &n)
[ac6e45]105{
106        this->fNormal=NULL;
[4199b3]107        this->interiorPoint=NULL;                       
[ac6e45]108        this->UCN=0;
109        this->codim2Ptr=NULL;
[ec6c52]110        if(n==2)
[ac6e45]111        {
112                this->codim=n;
[4199b3]113        }//NOTE Handle exception here!                 
[ac6e45]114        this->numCodim2Facets=0;
[5ece85]115        this->numRays=0;
[ac6e45]116        this->flipGB=NULL;
117        this->next=NULL;
118        this->prev=NULL;
119        this->flipRing=NULL;
120        this->isFlippable=FALSE;
121}
[4199b3]122               
[ac6e45]123/** \brief The copy constructor
[f5a3167]124* By default only copies the fNormal, f2Normals and UCN
[ac6e45]125*/
126facet::facet(const facet& f)
127{
[590f813]128        this->fNormal=iv64Copy(f.fNormal);
[f5a3167]129        this->UCN=f.UCN;
130        this->isFlippable=f.isFlippable;
[4199b3]131        //Needed for flip2
132        //NOTE ONLY REFERENCE
133        this->interiorPoint=iv64Copy(f.interiorPoint);//only referencing is prolly not the best thing to do in a copy constructor
[f5a3167]134        facet* f2Copy;
135        f2Copy=f.codim2Ptr;
136        facet* f2Act;
137        f2Act=this->codim2Ptr;
138        while(f2Copy!=NULL)
139        {
[f91fddc]140                if(f2Act==NULL
141#ifndef NDEBUG
142  #if SIZEOF_LONG==8
143                                 || f2Act==(facet*)0xfefefefefefefefe
144  #elif SIZEOF_LONG==4
145                                 || f2Act==(facet*)0xfefefefe
146  #endif
147#endif
148                  )
[f5a3167]149                {
150                        f2Act=new facet(2);
[4199b3]151                        this->codim2Ptr=f2Act;                 
[f5a3167]152                }
153                else
154                {
155                        facet* marker;
156                        marker = f2Act;
157                        f2Act->next = new facet(2);
158                        f2Act = f2Act->next;
159                        f2Act->prev = marker;
160                }
[590f813]161                int64vec *f2Normal;
[110e721]162                f2Normal = f2Copy->getFacetNormal();
[4199b3]163//              f2Act->setFacetNormal(f2Copy->getFacetNormal());
[110e721]164                f2Act->setFacetNormal(f2Normal);
165                delete f2Normal;
[f5a3167]166                f2Act->setUCN(f2Copy->getUCN());
167                f2Copy = f2Copy->next;
[4199b3]168        }       
[ac6e45]169}
[be7ab3f]170
171/** \brief Shallow copy constructor for facets
[73809b]172* We only need the interior point for equality testing
[be7ab3f]173*/
[73809b]174facet* facet::shallowCopy(const facet& f)
175{
176        facet *res = new facet();
[590f813]177        res->fNormal=(int64vec * const)f.fNormal;
[73809b]178        res->UCN=f.UCN;
179        res->isFlippable=f.isFlippable;
[590f813]180        res->interiorPoint=(int64vec * const)f.interiorPoint;
[73809b]181        res->codim2Ptr=(facet * const)f.codim2Ptr;
182        res->prev=NULL;
183        res->next=NULL;
184        res->flipGB=NULL;
[147167f]185        res->flipRing=NULL;
[73809b]186        return res;
187}
188
189void facet::shallowDelete()
190{
[96e408]191#ifndef NDEBUG
[4199b3]192//      printf("shallowdel@UCN %i\n", this->getUCN());
[147167f]193#endif
[73809b]194        this->fNormal=NULL;
[4199b3]195//      this->UCN=0;
[73809b]196        this->interiorPoint=NULL;
197        this->codim2Ptr=NULL;
198        this->prev=NULL;
199        this->next=NULL;
200        this->flipGB=NULL;
[147167f]201        this->flipRing=NULL;
[4199b3]202        assert(this->fNormal==NULL);   
203//      delete(this);
[73809b]204}
[4199b3]205               
[ac6e45]206/** The default destructor */
207facet::~facet()
208{
[96e408]209#ifndef NDEBUG
[4199b3]210//      printf("~facet@UCN %i\n",this->getUCN());
[147167f]211#endif
[1d37196]212        if(this->fNormal!=NULL)
213                delete this->fNormal;
214        if(this->interiorPoint!=NULL)
215                delete this->interiorPoint;
[963eeb]216        /* Cleanup the codim2-structure */
[4199b3]217//      if(this->codim==2)
218//      {
219//              facet *codim2Ptr;
220//              codim2Ptr = this->codim2Ptr;
221//              while(codim2Ptr!=NULL)
222//              {
223//                      if(codim2Ptr->fNormal!=NULL)
224//                      {
225//                              delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
226//                              codim2Ptr = codim2Ptr->next;
227//                      }
228//              }
229//      }
230        //The rays are stored in the cone!
[dcf8b4b]231        if(this->flipGB!=NULL)
232                idDelete((ideal *)&this->flipGB);
[4199b3]233//      if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
234//              rDelete(this->flipRing); //See vol II/134
235//      this->flipRing=NULL;
[dcf8b4b]236        this->prev=NULL;
237        this->next=NULL;
[ac6e45]238}
[4199b3]239               
[590f813]240inline const int64vec *facet::getRef2FacetNormal() const
[85583b]241{
242        return(this->fNormal);
[4199b3]243}       
[73809b]244
245/** Equality check for facets based on unique interior points*/
[ec6c52]246static bool areEqual2(facet* f, facet *g)
[be7ab3f]247{
248#ifdef gfanp
249        gcone::numberOfFacetChecks++;
250        timeval start, end;
251        gettimeofday(&start, 0);
252#endif
[73809b]253        bool res = TRUE;
[590f813]254        const int64vec *fIntP = f->getRef2InteriorPoint();
255        const int64vec *gIntP = g->getRef2InteriorPoint();
[1f637e]256        for(int ii=0;ii<(currRing->N);ii++)
[73809b]257        {
258                if( (*fIntP)[ii] != (*gIntP)[ii] )
259                {
260                        res=FALSE;
261                        break;
262                }
263        }
[4199b3]264//      if( fIntP->compare(gIntP)!=0) res=FALSE;
[be7ab3f]265#ifdef gfanp
266        gettimeofday(&end, 0);
[11a7dc]267        gcone::t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
[4199b3]268#endif 
[73809b]269        return res;
[be7ab3f]270}
271
[276932]272/** \brief Comparison of facets
273 * called from enqueueNewFacets
274* The facet normals are primitve vectors since we call gcone::normalize() on all cones.
275* Hence it should suffice to check whether facet normal f equals minus facet normal s.
276* If so we check the extremal rays
[85583b]277*
[590f813]278* BEWARE: It would be better to use const int64vec* but that will lead to call something like
279* int foo=((int64vec*)f2Normal)->compare((int64vec*)s2Normal) resulting in much higher memory usage
[276932]280*/
[ec6c52]281static bool areEqual(facet *f, facet *s)
[ac6e45]282{
[4c4a80]283#ifdef gfanp
284        gcone::numberOfFacetChecks++;
[85583b]285        timeval start, end;
286        gettimeofday(&start, 0);
[4c4a80]287#endif
[ac6e45]288        bool res = TRUE;
[f5a3167]289        int notParallelCtr=0;
290        int ctr=0;
[4199b3]291        const int64vec* fNormal; //No new since iv64Copy and therefore getFacetNormal return a new
[590f813]292        const int64vec* sNormal;
[e58dbb]293        fNormal = f->getRef2FacetNormal();
294        sNormal = s->getRef2FacetNormal();
[76cfef]295#include <misc/intvec.h>
[4199b3]296        //Do not need parallelity. Too time consuming
297//      if(!isParallel(*fNormal,*sNormal))
298//      if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
299 //             notParallelCtr++;
300//      else//parallelity, so we check the codim2-facets
[590f813]301        int64vec *fNRef=const_cast<int64vec*>(fNormal);
302        int64vec *sNRef=const_cast<int64vec*>(sNormal);
[5e4025]303        if(isParallel(*fNRef,*sNRef))
[4199b3]304//      if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
[ac6e45]305        {
[f5a3167]306                facet* f2Act;
307                facet* s2Act;
[4199b3]308                f2Act = f->codim2Ptr;           
[f5a3167]309                ctr=0;
310                while(f2Act!=NULL)
[008e81]311                {
[590f813]312                        const int64vec* f2Normal;
[e58dbb]313                        f2Normal = f2Act->getRef2FacetNormal();
[4199b3]314//                      int64vec *f2Ref=const_cast<int64vec*>(f2Normal);
[f5a3167]315                        s2Act = s->codim2Ptr;
316                        while(s2Act!=NULL)
[008e81]317                        {
[590f813]318                                const int64vec* s2Normal;
[e58dbb]319                                s2Normal = s2Act->getRef2FacetNormal();
[4199b3]320//                              bool foo=areEqual(f2Normal,s2Normal);
321//                              int64vec *s2Ref=const_cast<int64vec*>(s2Normal);
[4c4a80]322                                int foo=f2Normal->compare(s2Normal);
323                                if(foo==0)
[f5a3167]324                                        ctr++;
[dcf8b4b]325                                s2Act = s2Act->next;
[f5a3167]326                        }
327                        f2Act = f2Act->next;
[4199b3]328                }               
[ac6e45]329        }
[f5a3167]330        if(ctr==f->numCodim2Facets)
331                res=TRUE;
332        else
[4c4a80]333        {
334#ifdef gfanp
335                gcone::parallelButNotEqual++;
336#endif
337                res=FALSE;
338        }
[85583b]339#ifdef gfanp
340        gettimeofday(&end, 0);
[11a7dc]341        gcone::t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
[85583b]342#endif
[ac6e45]343        return res;
[4199b3]344}       
345               
[590f813]346/** Stores the facet normal \param int64vec*/
347inline void facet::setFacetNormal(int64vec *iv)
[ac6e45]348{
[2f72c55]349        if(this->fNormal!=NULL)
350                delete this->fNormal;
[4199b3]351        this->fNormal = iv64Copy(iv);                   
[ac6e45]352}
[4199b3]353               
354/** Hopefully returns the facet normal
[590f813]355* Mind: iv64Copy returns a new int64vec, so use this in the following way:
356* int64vec *iv;
[72e467]357* iv = this->getFacetNormal();
358* [...]
359* delete(iv);
360*/
[590f813]361inline int64vec *facet::getFacetNormal() const
[4199b3]362{                               
[590f813]363        return iv64Copy(this->fNormal);
[ac6e45]364}
[73809b]365
[ac6e45]366/** Method to print the facet normal*/
[ce41a2e]367inline void facet::printNormal() const
[ac6e45]368{
369        fNormal->show();
370}
[4199b3]371               
[ac6e45]372/** Store the flipped GB*/
[a23a297]373inline void facet::setFlipGB(ideal I)
[ac6e45]374{
375        this->flipGB=idCopy(I);
376}
[4199b3]377               
[4c4a80]378/** Returns a pointer to the flipped GB
[d4f1b95]379Seems not be used
380Anyhow idCopy would make sense here.
381*/
[a23a297]382inline ideal facet::getFlipGB()
[ac6e45]383{
384        return this->flipGB;
385}
[4199b3]386               
[ac6e45]387/** Print the flipped GB*/
[a23a297]388inline void facet::printFlipGB()
[ac6e45]389{
[3b9cf6]390#ifndef NDEBUG
[ac6e45]391        idShow(this->flipGB);
[3b9cf6]392#endif
[ac6e45]393}
[4199b3]394               
[ac6e45]395/** Set the UCN */
[a23a297]396inline void facet::setUCN(int n)
[ac6e45]397{
398        this->UCN=n;
399}
[4199b3]400               
401/** \brief Get the UCN
[ac6e45]402* Returns the UCN iff this != NULL, else -1
403*/
[a23a297]404inline int facet::getUCN()
[ac6e45]405{
[1f3450]406#ifndef NDEBUG
407  #if SIZEOF_LONG==8
408        if((this!=NULL && this!=(facet * const)0xfbfbfbfbfbfbfbfb))
[c82d55]409  #elif SIZEOF_LONG==4
[1f3450]410        if((this!=NULL && this!=(facet * const)0xfbfbfbfb))
411  #endif
412#endif
413#ifdef NDEBUG
414        if(this!=NULL)
[4199b3]415#endif                 
[ac6e45]416                return this->UCN;
417        else
418                return -1;
419}
[4199b3]420               
[ac6e45]421/** Store an interior point of the facet */
[590f813]422inline void facet::setInteriorPoint(int64vec *iv)
[ac6e45]423{
[2f72c55]424        if(this->interiorPoint!=NULL)
425                delete this->interiorPoint;
[590f813]426        this->interiorPoint = iv64Copy(iv);
[ac6e45]427}
[4199b3]428               
[4c4a80]429/** Returns a pointer to this->interiorPoint
[590f813]430* MIND: iv64Copy returns a new int64vec
[4c4a80]431* @see facet::getFacetNormal
432*/
[590f813]433inline int64vec *facet::getInteriorPoint()
[ac6e45]434{
[590f813]435        return iv64Copy(this->interiorPoint);
[73809b]436}
437
[590f813]438inline const int64vec *facet::getRef2InteriorPoint()
[73809b]439{
440        return (this->interiorPoint);
441}
[4199b3]442               
[ac6e45]443/** \brief Debugging function
444* prints the facet normal an all (codim-2)-facets that belong to it
445*/
[be7ab3f]446volatile void facet::fDebugPrint()
[e58dbb]447{
448  #ifndef NDEBUG
[4199b3]449        facet *codim2Act;                       
[ac6e45]450        codim2Act = this->codim2Ptr;
[590f813]451        int64vec *fNormal;
[4199b3]452        fNormal = this->getFacetNormal();       
[11a7dc]453        printf("=======================\n");
454        printf("Facet normal = (");fNormal->show(1,1);printf(")\n");
455        printf("-----------------------\n");
456        printf("Codim2 facets:\n");
[ac6e45]457        while(codim2Act!=NULL)
458        {
[590f813]459                int64vec *f2Normal;
[ac6e45]460                f2Normal = codim2Act->getFacetNormal();
[11a7dc]461                printf("(");f2Normal->show(1,0);printf(")\n");
[ac6e45]462                codim2Act = codim2Act->next;
[dcf8b4b]463                delete f2Normal;
[ac6e45]464        }
[11a7dc]465        printf("=======================\n");
[e58dbb]466        delete fNormal;
467  #endif
468}
[8bdaab]469
[d3ce3f1]470
[fe919c]471/**
[90adba8]472*\brief Implements the cone structure
473*
474* A cone is represented by a linked list of facet normals
475* @see facet
[fe919c]476*/
[1d9469]477
478
[4199b3]479/** \brief Default constructor.
[1d9469]480 *
481 * Initialises this->next=NULL and this->facetPtr=NULL
482 */
483gcone::gcone()
[8bdaab]484{
[1d9469]485        this->next=NULL;
[c5940e]486        this->prev=NULL;
[4199b3]487        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
[1d9469]488        this->baseRing=currRing;
489        this->counter++;
[4199b3]490        this->UCN=this->counter;                       
[1d9469]491        this->numFacets=0;
492        this->ivIntPt=NULL;
[68eb0c]493        this->gcRays=NULL;
[1d9469]494}
[4199b3]495               
[1d9469]496/** \brief Constructor with ring and ideal
497 *
[4199b3]498 * This constructor takes the root ring and the root ideal as parameters and stores
[1d9469]499 * them in the private members gcone::rootRing and gcone::inputIdeal
[e2e3ad]500 * This constructor is only called once in the computation of the Gröbner fan,
501 * namely for the very first cone. Therefore pred is set to 1.
502 * Might set it to this->UCN though...
[1d9469]503 * Since knowledge of the root ring is only needed when using reverse search,
504 * this constructor is not needed when using the "second" method
505*/
506gcone::gcone(ring r, ideal I)
507{
508        this->next=NULL;
[c5940e]509        this->prev=NULL;
[4199b3]510        this->facetPtr=NULL;                   
[1d9469]511        this->inputIdeal=I;
512        this->baseRing=currRing;
513        this->counter++;
[e2e3ad]514        this->UCN=this->counter;
515        this->pred=1;
[1d9469]516        this->numFacets=0;
517        this->ivIntPt=NULL;
[68eb0c]518        this->gcRays=NULL;
[1d9469]519}
[4199b3]520               
521/** \brief Copy constructor
[1d9469]522 *
523 * Copies a cone, sets this->gcBasis to the flipped GB
524 * Call this only after a successful call to gcone::flip which sets facet::flipGB
[4199b3]525*/             
[1d9469]526gcone::gcone(const gcone& gc, const facet &f)
527{
528        this->next=NULL;
[4199b3]529//      this->prev=(gcone *)&gc; //comment in to get a tree
[c5940e]530        this->prev=NULL;
[e58dbb]531        this->numVars=gc.numVars;
[1d9469]532        this->counter++;
[e2e3ad]533        this->UCN=this->counter;
[a23a297]534        this->pred=gc.UCN;
[1d9469]535        this->facetPtr=NULL;
[6123d46]536        this->gcBasis=idrCopyR(f.flipGB, f.flipRing);
[4199b3]537//      this->inputIdeal=idCopy(this->gcBasis);
[1d9469]538        this->baseRing=rCopy(f.flipRing);
539        this->numFacets=0;
540        this->ivIntPt=NULL;
[68eb0c]541        this->gcRays=NULL;
[1d9469]542}
[4199b3]543               
544/** \brief Default destructor
[1d9469]545*/
546gcone::~gcone()
547{
[73809b]548#ifndef NDEBUG
549  #if SIZEOF_LONG==8
[147167f]550        if( ( this->gcBasis!=(ideal)(0xfbfbfbfbfbfbfbfb) ) && (this->gcBasis!=NULL) )
[73809b]551                idDelete((ideal*)&this->gcBasis);
552  #elif SIZEOF_LONG!=8
553        if(this->gcBasis!=(ideal)0xfbfbfbfb)
554                idDelete((ideal *)&this->gcBasis);
555  #endif
556#else
[bec3b0]557        if(this->gcBasis!=NULL)
558                idDelete((ideal *)&this->gcBasis);
[73809b]559#endif
[4199b3]560//      idDelete((ideal *)&this->gcBasis);
561//      if(this->inputIdeal!=NULL)
562//              idDelete((ideal *)&this->inputIdeal);
563//      if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
564//              rDelete(this->rootRing);
[1f3450]565        if(this->UCN!=1 && this->baseRing!=NULL)
566                rDelete(this->baseRing);
[c5940e]567        facet *fAct;
568        facet *fDel;
569        /*Delete the facet structure*/
570        fAct=this->facetPtr;
571        fDel=fAct;
572        while(fAct!=NULL)
573        {
574                fDel=fAct;
575                fAct=fAct->next;
576                delete fDel;
577        }
[f07b38c]578        this->counter--;
[4199b3]579        //should be deleted in noRevS
580//      dd_FreeMatrix(this->ddFacets);
581        //dd_FreeMatrix(this->ddFacets);
[68eb0c]582        for(int ii=0;ii<this->numRays;ii++)
583                delete(gcRays[ii]);
584        omFree(gcRays);
[4199b3]585}                       
[73809b]586
587/** Returns the number of cones existing at the time*/
[a23a297]588inline int gcone::getCounter()
[0d278cd]589{
590        return this->counter;
591}
[4199b3]592       
[1d9469]593/** \brief Set the interior point of a cone */
[590f813]594inline void gcone::setIntPoint(int64vec *iv)
[1d9469]595{
[2f72c55]596        if(this->ivIntPt!=NULL)
597                delete this->ivIntPt;
[590f813]598        this->ivIntPt=iv64Copy(iv);
[1d9469]599}
[4199b3]600               
[73809b]601/** \brief Returns either a physical copy the interior point of a cone or just a reference to it.*/
[590f813]602inline int64vec *gcone::getIntPoint(bool shallow)
[1d9469]603{
[73809b]604        if(shallow==TRUE)
605                return this->ivIntPt;
606        else
[590f813]607                return iv64Copy(this->ivIntPt);
[1d9469]608}
[4199b3]609               
[1d9469]610/** \brief Print the interior point */
[a23a297]611inline void gcone::showIntPoint()
[1d9469]612{
613        ivIntPt->show();
614}
[4199b3]615               
[1d9469]616/** \brief Print facets
617 * This is mainly for debugging purposes. Usually called from within gdb
618 */
[5ece85]619volatile void gcone::showFacets(const short codim)
[1d9469]620{
[e58dbb]621  #ifndef NDEBUG
[5ece85]622        facet *f=this->facetPtr;
623        facet *f2=NULL;
624        if(codim==2)
625                f2=this->facetPtr->codim2Ptr;
[1d9469]626        while(f!=NULL)
627        {
[590f813]628                int64vec *iv;
[1d9469]629                iv = f->getFacetNormal();
[11a7dc]630                printf("(");iv->show(1,0);
[1d9469]631                if(f->isFlippable==FALSE)
[11a7dc]632                        printf(")* ");
[1d9469]633                else
[11a7dc]634                        printf(") ");
[dcf8b4b]635                delete iv;
[5ece85]636                if(codim==2)
637                {
638                        f2=f->codim2Ptr;
639                        while(f2!=NULL)
640                        {
[11a7dc]641                                printf("[");f2->getFacetNormal()->show(1,0);printf("]");
[5ece85]642                                f2 = f2->next;
643                        }
[11a7dc]644                        printf("\n");
[5ece85]645                }
[1d9469]646                f=f->next;
647        }
[11a7dc]648        printf("\n");
[e58dbb]649  #endif
[1d9469]650}
[4199b3]651               
[1d9469]652/** For debugging purposes only */
[ec6c52]653static volatile void showSLA(facet &f)
[1d9469]654{
[e58dbb]655  #ifndef NDEBUG
[1d9469]656        facet *fAct;
657        fAct = &f;
[1f3450]658        if(fAct!=NULL)
[1d9469]659        {
[1f3450]660                facet *codim2Act;
[1d9469]661                codim2Act = fAct->codim2Ptr;
[4199b3]662               
[11a7dc]663                printf("\n");
[1f3450]664                while(fAct!=NULL)
[493699e]665                {
[4199b3]666                        int64vec *fNormal;             
[1f3450]667                        fNormal=fAct->getFacetNormal();
[11a7dc]668                        printf("(");fNormal->show(1,0);
[1f3450]669                        if(fAct->isFlippable==TRUE)
[11a7dc]670                                printf(") ");
[1f3450]671                        else
[11a7dc]672                                printf(")* ");
[1f3450]673                        delete fNormal;
674                        codim2Act = fAct->codim2Ptr;
[11a7dc]675                        printf(" Codim2: ");
[1f3450]676                        while(codim2Act!=NULL)
677                        {
[590f813]678                                int64vec *f2Normal;
[1f3450]679                                f2Normal = codim2Act->getFacetNormal();
[11a7dc]680                                printf("(");f2Normal->show(1,0);printf(") ");
[1f3450]681                                delete f2Normal;
682                                codim2Act = codim2Act->next;
683                        }
[63a2f8]684                        printf("UCN = %i\n",fAct->getUCN());
[1f3450]685                        fAct = fAct->next;
[493699e]686                }
[1d9469]687        }
[e58dbb]688  #endif
[1d9469]689}
[4199b3]690               
[ec6c52]691static void idDebugPrint(const ideal &I)
[1d9469]692{
[e58dbb]693  #ifndef NDEBUG
[1d9469]694        int numElts=IDELEMS(I);
[11a7dc]695        printf("Ideal with %i generators\n", numElts);
696        printf("Leading terms: ");
[1d9469]697        for (int ii=0;ii<numElts;ii++)
698        {
699                pWrite0(pHead(I->m[ii]));
[11a7dc]700                printf(",");
[1d9469]701        }
[11a7dc]702        printf("\n");
[e58dbb]703  #endif
[1d9469]704}
[913422]705
[ec6c52]706static void invPrint(const ideal &I)
[7c0ec6]707{
[4199b3]708//      int numElts=IDELEMS(I);
709//      cout << "inv = ";
710//      for(int ii=0;ii<numElts;ii++);
711//      {
712//              pWrite0(pHead(I->m[ii]));
713//              cout << ",";
714//      }
715//      cout << endl;
[7c0ec6]716}
717
[ec6c52]718static bool isMonomial(const ideal &I)
[913422]719{
720        bool res = TRUE;
721        for(int ii=0;ii<IDELEMS(I);ii++)
722        {
723                if(pLength((poly)I->m[ii])>1)
724                {
725                        res = FALSE;
726                        break;
[276932]727                }
[913422]728        }
729        return res;
730}
[4199b3]731               
[1d9469]732/** \brief Set gcone::numFacets */
[a23a297]733inline void gcone::setNumFacets()
[1d9469]734{
735}
[4199b3]736               
[1d9469]737/** \brief Get gcone::numFacets */
[a23a297]738inline int gcone::getNumFacets()
[1d9469]739{
740        return this->numFacets;
741}
[4199b3]742               
[a23a297]743inline int gcone::getUCN()
[1d9469]744{
[4199b3]745        if( this!=NULL)// && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
[1d9469]746                return this->UCN;
747        else
748                return -1;
749}
[e2e3ad]750
[a23a297]751inline int gcone::getPredUCN()
[e2e3ad]752{
753        return this->pred;
754}
[73809b]755/** Returns a copy of the this->baseRing */
756inline ring gcone::getBaseRing()
757{
758        return rCopy(this->baseRing);
759}
[e2e3ad]760
[63a2f8]761inline void gcone::setBaseRing(ring r)
762{
763        this->baseRing=rCopy(r);
764}
765
[147167f]766inline ring gcone::getRef2BaseRing()
767{
768        return this->baseRing;
769}
770
[1d9469]771/** \brief Compute the normals of the cone
772 *
773 * This method computes a representation of the cone in terms of facet normals. It takes an ideal
774 * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
775 * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
[4199b3]776 * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
[1d9469]777 * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
778 * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
779 * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
780 *
781 * Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
782 * an interior point of the cone.
783                 */
[ec6c52]784void gcone::getConeNormals(const ideal &I, bool compIntPoint)
[1d9469]785{
[42660f]786#ifdef gfanp
787        timeval start, end;
788        gettimeofday(&start, 0);
789#endif
[1d9469]790        poly aktpoly;
[4199b3]791        int rows;                       // will contain the dimensions of the ineq matrix - deprecated by
[1d9469]792        dd_rowrange ddrows;
793        dd_colrange ddcols;
[4199b3]794        dd_rowset ddredrows;            // # of redundant rows in ddineq
795        dd_rowset ddlinset;             // the opposite
796        dd_rowindex ddnewpos=NULL;                // all to make dd_Canonicalize happy
797        dd_NumberType ddnumb=dd_Integer; //Number type
798        dd_ErrorType dderr=dd_NoError;         
799        //Compute the # inequalities i.e. rows of the matrix
800        rows=0; //Initialization
[1d9469]801        for (int ii=0;ii<IDELEMS(I);ii++)
802        {
[4199b3]803//              aktpoly=(poly)I->m[ii];
804//              rows=rows+pLength(aktpoly)-1;
[73809b]805                rows=rows+pLength((poly)I->m[ii])-1;
[1d9469]806        }
[e2e3ad]807
[4199b3]808        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
[1d9469]809        ddrows=rows;
[e2e3ad]810        ddcols=this->numVars;
[4199b3]811        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
812        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
813               
814        // We loop through each g\in GB and compute the resulting inequalities
[1d9469]815        for (int i=0; i<IDELEMS(I); i++)
816        {
[4199b3]817                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
818                //simpler version of storing expvect diffs
[e2e3ad]819                int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
[0d278cd]820                pGetExpV(aktpoly,leadexpv);
[63a2f8]821                poly pNextTerm=aktpoly;
822                while(pNext(pNextTerm)/*pNext(aktpoly)*/!=NULL)
[1d9469]823                {
[63a2f8]824                        pNextTerm/*aktpoly*/=pNext(pNextTerm);
[42660f]825                        int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
[4199b3]826                        pGetExpV(pNextTerm,tailexpv);                   
[e2e3ad]827                        for(int kk=1;kk<=this->numVars;kk++)
[4199b3]828                        {                               
[0d278cd]829                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
[1d9469]830                        }
[0d278cd]831                        aktmatrixrow += 1;
[42660f]832                        omFree(tailexpv);
[4199b3]833                }               
834                omFree(leadexpv);       
835        } //for
[b45f583]836#if true
[7c0ec6]837        /*Let's make the preprocessing here. This could already be done in the above for-loop,
838        * but for a start it is more convenient here.
839        * We check the necessary condition of FJT p.18
840        * Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
841        */
[4199b3]842//      ideal initialForm=idInit(IDELEMS(I),1);
843//      int64vec *gamma=new int64vec(this->numVars);
[7c0ec6]844        int falseGammaCounter=0;
845        int *redRowsArray=NULL;
846        int num_alloc=0;
[4199b3]847        int num_elts=0; 
[7c0ec6]848        for(int ii=0;ii<ddineq->rowsize;ii++)
849        {
[f91fddc]850                ideal initialForm=idInit(IDELEMS(I),I->rank);
[4199b3]851                //read row ii into gamma
852//              int64 tmp;
[590f813]853                int64vec *gamma=new int64vec(this->numVars);
[7c0ec6]854                for(int jj=1;jj<=this->numVars;jj++)
855                {
[f91fddc]856                        int64 tmp;
[590f813]857                        tmp=(int64)mpq_get_d(ddineq->matrix[ii][jj]);
858                        (*gamma)[jj-1]=(int64)tmp;
[7c0ec6]859                }
860                computeInv((ideal&)I,initialForm,*gamma);
[dd6b1d]861                delete gamma;
[4199b3]862                //Create leading ideal
[7c0ec6]863                ideal L=idInit(IDELEMS(initialForm),1);
864                for(int jj=0;jj<IDELEMS(initialForm);jj++)
865                {
[f91fddc]866                        poly p=pHead(initialForm->m[jj]);
867                        L->m[jj]=pCopy(/*pHead(initialForm->m[jj]))*/p);
868                        pDelete(&p);
[4199b3]869                }               
870               
[871f01]871                LObject *P = new sLObject();//TODO What's the difference between sLObject and LObject?
[7c0ec6]872                memset(P,0,sizeof(LObject));
[b45f583]873
[7c0ec6]874                for(int jj=0;jj<=IDELEMS(initialForm)-2;jj++)
875                {
876                        bool isMaybeFacet=FALSE;
[4199b3]877                        P->p1=initialForm->m[jj];       //build spolys of initialForm in_v
[b45f583]878
[7c0ec6]879                        for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
880                        {
881                                P->p2=initialForm->m[kk];
[4199b3]882                                ksCreateSpoly(P);                                                       
883                                if(P->p!=NULL)  //spoly non zero=?
884                                {       
[871f01]885                                        poly p;//NOTE Don't use pInit here. Evil memleak will follow
[dd6b1d]886                                        poly q;
887                                        poly pDel,qDel;
[7c0ec6]888                                        p=pCopy(P->p);
[4199b3]889                                        q=pHead(p);     //Monomial q
[dd6b1d]890                                        pDelete(&q);
891                                        pDel=p; qDel=q;
[7c0ec6]892                                        isMaybeFacet=FALSE;
[4199b3]893                                        //TODO: Suffices to check LTs here
[7c0ec6]894                                        while(p!=NULL)
[4199b3]895                                        {                                               
[63a2f8]896                                                q=pHead(p);
[7c0ec6]897                                                for(int ll=0;ll<IDELEMS(L);ll++)
898                                                {
899                                                        if(pLmEqual(L->m[ll],q) || pDivisibleBy(L->m[ll],q))
[4199b3]900                                                        {                                                       
[7c0ec6]901                                                                isMaybeFacet=TRUE;
[4199b3]902                                                                break;//for
[7c0ec6]903                                                        }
904                                                }
[f91fddc]905                                                pDelete(&q);
[7c0ec6]906                                                if(isMaybeFacet==TRUE)
[f91fddc]907                                                {
[4199b3]908                                                        break;//while(p!=NULL)
[7c0ec6]909                                                }
[dd6b1d]910                                                p=pNext(p);
[4199b3]911                                        }//while
912//                                      pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
[f91fddc]913                                        if(q!=NULL) pDelete(&q);
[dd6b1d]914                                        pDelete(&pDel);
915                                        pDelete(&qDel);
[7c0ec6]916                                        if(isMaybeFacet==FALSE)
917                                        {
[4199b3]918                                                dd_set_si(ddineq->matrix[ii][0],1);                                             
919//                                              if(num_alloc==0)
920//                                                      num_alloc += 1;
921//                                              else                                           
922//                                                      num_alloc += 1;
[f91fddc]923                                                if(num_alloc==num_elts) num_alloc==0 ? num_alloc=1 : num_alloc*=2;
[4199b3]924                                               
[1d37196]925                                                void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
926                                                if(!tmp)
927                                                {
[b45f583]928                                                        WerrorS("Woah dude! Couldn't realloc memory\n");
[1d37196]929                                                        exit(-1);
930                                                }
931                                                redRowsArray = (int*)tmp;
[7c0ec6]932                                                redRowsArray[num_elts]=ii;
933                                                num_elts++;
[4199b3]934                                                //break;//for(int kk, since we have found one that is not in L 
935                                                goto _start;    //mea culpa, mea culpa, mea maxima culpa
[7c0ec6]936                                        }
[4199b3]937                                }//if(P->p!=NULL)
[dd6b1d]938                                pDelete(&(P->p));
[4199b3]939                        }//for k
940                }//for jj
[b45f583]941                _start:;
[dcf8b4b]942                idDelete(&L);
[7c0ec6]943                delete P;
[dcf8b4b]944                idDelete(&initialForm);
[4199b3]945        }//for(ii<ddineq-rowsize
946//      delete gamma;
947        int offset=0;//needed for correction of redRowsArray[ii]
[96e408]948#ifndef NDEBUG
[f91fddc]949        printf("Removed %i of %i in preprocessing step\n",num_elts,ddineq->rowsize);
950#endif
[7c0ec6]951        for( int ii=0;ii<num_elts;ii++ )
952        {
[4199b3]953                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
[7c0ec6]954                offset++;
[4199b3]955        }       
956        free(redRowsArray);//NOTE May crash on some machines.
[2e06300]957        /*And now for the strictly positive rows
958        * Doesn't gain significant speedup
959        */
960        /*int *posRowsArray=NULL;
961        num_alloc=0;
962        num_elts=0;
963        for(int ii=0;ii<ddineq->rowsize;ii++)
964        {
[590f813]965                int64vec *ivPos = new int64vec(this->numVars);
[2e06300]966                for(int jj=0;jj<this->numVars;jj++)
967                        (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
968                bool isStrictlyPos=FALSE;
[4199b3]969                int posCtr=0;           
[2e06300]970                for(int jj=0;jj<this->numVars;jj++)
971                {
[590f813]972                        int64vec *ivCanonical = new int64vec(this->numVars);
[2e06300]973                        jj==0 ? (*ivCanonical)[ivPos->length()-1]=1 : (*ivCanonical)[jj-1]=1;
974                        if(dotProduct(*ivCanonical,*ivPos)!=0)
975                        {
976                                if ((*ivPos)[jj]>=0)
[4199b3]977                                {                               
978                                        posCtr++;                               
[2e06300]979                                }
[4199b3]980                        }                       
[2e06300]981                        delete ivCanonical;
982                }
983                if(posCtr==ivPos->length())
984                        isStrictlyPos=TRUE;
985                if(isStrictlyPos==TRUE)
986                {
987                        if(num_alloc==0)
988                                num_alloc += 1;
989                        else
990                                num_alloc += 1;
991                        void *tmp = realloc(posRowsArray,(num_alloc*sizeof(int)));
992                        if(!tmp)
993                        {
994                                WerrorS("Woah dude! Couldn't realloc memory\n");
995                                exit(-1);
996                        }
997                        posRowsArray = (int*)tmp;
998                        posRowsArray[num_elts]=ii;
[4199b3]999                        num_elts++;     
[2e06300]1000                }
1001                delete ivPos;
1002        }
1003        offset=0;
1004        for(int ii=0;ii<num_elts;ii++)
1005        {
1006                dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
1007                offset++;
1008        }
1009        free(posRowsArray);*/
[7c0ec6]1010#endif
[c5d8dd]1011
[4199b3]1012        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);       
1013        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
[1d9469]1014        ddcols = ddineq->colsize;
[4199b3]1015       
[4c4a80]1016        this->ddFacets = dd_CopyMatrix(ddineq);
[4199b3]1017                       
1018        /*Write the normals into class facet*/ 
1019        facet *fAct;    //pointer to active facet       
[1d9469]1020        int numNonFlip=0;
1021        for (int kk = 0; kk<ddrows; kk++)
1022        {
[4199b3]1023                int64 ggT=1;//NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
1024                int64vec *load = new int64vec(this->numVars);//int64vec to store a single facet normal that will then be stored via setFacetNormal
[1d9469]1025                for (int jj = 1; jj <ddcols; jj++)
1026                {
[f91fddc]1027                        int64 val;
1028                        val = (int64)mpq_get_d(ddineq->matrix[kk][jj]);
[4199b3]1029                        (*load)[jj-1] = val;    //store typecasted entry at pos jj-1 of load
[f91fddc]1030                        ggT = int64gcd(ggT,/*(int64&)foo*/val);
[4199b3]1031                }//for (int jj = 1; jj <ddcols; jj++)
[73809b]1032                if(ggT>1)
1033                {
1034                        for(int ll=0;ll<this->numVars;ll++)
[4199b3]1035                                (*load)[ll] /= ggT;//make primitive vector
[73809b]1036                }
[4c4a80]1037                /*Quick'n'dirty hack for flippability. Executed only if gcone::hasHomInput==FALSE
1038                * Otherwise every facet intersects the positive orthant
[4199b3]1039                */     
[4c4a80]1040                if(gcone::hasHomInput==FALSE)
1041                {
[4199b3]1042                        //TODO: No dP needed
[72e467]1043                        bool isFlip=FALSE;
[4c4a80]1044                        for(int jj = 0; jj<load->length(); jj++)
[1d9469]1045                        {
[4199b3]1046//                              int64vec *ivCanonical = new int64vec(load->length());
1047//                              (*ivCanonical)[jj]=1;
1048//                              if (dotProduct(*load,*ivCanonical)<0)   
1049//                              {
1050//                                      isFlip=TRUE;
1051//                                      break;  //URGHS
1052//                              }
1053//                              delete ivCanonical;
[85583b]1054                                if((*load)[jj]<0)
[72e467]1055                                {
1056                                        isFlip=TRUE;
[85583b]1057                                        break;
[4199b3]1058                                }                               
[72e467]1059                        }/*End of check for flippability*/
[4199b3]1060//                      if(iv64isStrictlyPositive(load))
1061//                              isFlip=TRUE;
[4c4a80]1062                        if(isFlip==FALSE)
[1d9469]1063                        {
[72e467]1064                                this->numFacets++;
1065                                numNonFlip++;
1066                                if(this->numFacets==1)
1067                                {
1068                                        facet *fRoot = new facet();
1069                                        this->facetPtr = fRoot;
[4199b3]1070                                        fAct = fRoot;                                                   
[72e467]1071                                }
1072                                else
1073                                {
1074                                        fAct->next = new facet();
1075                                        fAct = fAct->next;
1076                                }
1077                                fAct->isFlippable=FALSE;
1078                                fAct->setFacetNormal(load);
1079                                fAct->setUCN(this->getUCN());
[96e408]1080#ifndef NDEBUG
[4199b3]1081                                printf("Marking facet (");load->show(1,0);printf(") as non flippable\n");               
[b7510d]1082#endif
[72e467]1083                        }
[4c4a80]1084                        else
1085                        {
1086                                this->numFacets++;
1087                                if(this->numFacets==1)
1088                                {
1089                                        facet *fRoot = new facet();
1090                                        this->facetPtr = fRoot;
1091                                        fAct = fRoot;
1092                                }
1093                                else
1094                                {
1095                                        fAct->next = new facet();
1096                                        fAct = fAct->next;
1097                                }
1098                                fAct->isFlippable=TRUE;
1099                                fAct->setFacetNormal(load);
[4199b3]1100                                fAct->setUCN(this->getUCN());                                   
[4c4a80]1101                        }
[4199b3]1102                }//hasHomInput==FALSE
1103                else    //Every facet is flippable
1104                {       /*Now load should be full and we can call setFacetNormal*/                                     
[1d9469]1105                        this->numFacets++;
1106                        if(this->numFacets==1)
1107                        {
1108                                facet *fRoot = new facet();
1109                                this->facetPtr = fRoot;
1110                                fAct = fRoot;
1111                        }
1112                        else
1113                        {
1114                                fAct->next = new facet();
1115                                fAct = fAct->next;
1116                        }
1117                        fAct->isFlippable=TRUE;
1118                        fAct->setFacetNormal(load);
[4199b3]1119                        fAct->setUCN(this->getUCN());                                   
1120                }//if (isFlippable==FALSE)
1121                delete load;                           
1122        }//for (int kk = 0; kk<ddrows; kk++)
1123                       
1124        //In cases like I=<x-1,y-1> there are only non-flippable facets...
[1d9469]1125        if(numNonFlip==this->numFacets)
[e58dbb]1126        {
[e2e3ad]1127                WerrorS ("Only non-flippable facets. Terminating...\n");
[4199b3]1128//              exit(-1);//Bit harsh maybe...
[1d9469]1129        }
[4199b3]1130                       
[963eeb]1131        /*
[1d9469]1132        Now we should have a linked list containing the facet normals of those facets that are
1133        -irredundant
1134        -flipable
1135        Adressing is done via *facetPtr
[4199b3]1136        */                     
[1d9469]1137        if (compIntPoint==TRUE)
1138        {
[590f813]1139                int64vec *iv = new int64vec(this->numVars);
[1d9469]1140                dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1141                int jj=1;
1142                for (int ii=0;ii<=this->numVars;ii++)
1143                {
1144                        dd_set_si(posRestr->matrix[ii][jj],1);
[4199b3]1145                        jj++;                                                   
[1d9469]1146                }
1147                dd_MatrixAppendTo(&ddineq,posRestr);
[4199b3]1148                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
1149                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
[e2e3ad]1150                delete iv;
[dcf8b4b]1151                dd_FreeMatrix(posRestr);
[4199b3]1152        }       
1153        //Clean up but don't delete the return value!   
1154        //dd_FreeMatrix(ddineq);       
1155        set_free(ddredrows);//check
1156        set_free(ddlinset);//check
1157        //free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
[42660f]1158#ifdef gfanp
1159        gettimeofday(&end, 0);
1160        time_getConeNormals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1161#endif
[90adba8]1162
[4199b3]1163}//gcone::getConeNormals(ideal I)
1164               
[1d9469]1165/** \brief Compute the (codim-2)-facets of a given cone
1166 * This method is used during noRevS
1167 * Additionally we check whether the codim2-facet normal is strictly positive. Otherwise
1168 * the facet is marked as non-flippable.
1169 */
[ec6c52]1170void gcone::getCodim2Normals(const gcone &gc)
[1d9469]1171{
[42660f]1172#ifdef gfanp
1173        timeval start, end;
1174        gettimeofday(&start, 0);
1175#endif
[4199b3]1176        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
[1d9469]1177        facet *fAct;
[4199b3]1178        fAct = this->facetPtr;         
[1d9469]1179        facet *codim2Act;
[4199b3]1180        //codim2Act = this->facetPtr->codim2Ptr;
1181        dd_MatrixPtr ddineq;//,P,ddakt;
[1d9469]1182        dd_ErrorType err;
[4199b3]1183        //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
[1d9469]1184        ddineq = dd_CopyMatrix(gc.ddFacets);
1185        /*Now set appropriate linearity*/
[4199b3]1186        for (int ii=0; ii<this->numFacets; ii++)                       
1187        {       
[42660f]1188                dd_rowset impl_linset, redset;
1189                dd_rowindex newpos;
1190                dd_MatrixPtr ddakt;
1191                ddakt = dd_CopyMatrix(ddineq);
[4199b3]1192//              ddakt->representation=dd_Inequality;    //Not using this makes it faster. But why does the quick check below still work?
1193//              ddakt->representation=dd_Generator;
[42660f]1194                set_addelem(ddakt->linset,ii+1);/*Now set appropriate linearity*/
1195#ifdef gfanp
1196                timeval t_ddMC_start, t_ddMC_end;
1197                gettimeofday(&t_ddMC_start,0);
[4199b3]1198#endif                         
1199                //dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
[42660f]1200                dd_PolyhedraPtr ddpolyh;
[1d9469]1201                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
[4199b3]1202//              ddpolyh=dd_DDMatrix2Poly2(ddakt, dd_MaxCutoff, &err);
[42660f]1203                dd_MatrixPtr P;
[4199b3]1204                P=dd_CopyGenerators(ddpolyh);           
[42660f]1205                dd_FreePolyhedra(ddpolyh);
[4199b3]1206                //TODO Call for one cone , normalize - check equalities - plus lineality -done
[42660f]1207#ifdef gfanp
1208                gettimeofday(&t_ddMC_end,0);
1209                t_ddMC += (t_ddMC_end.tv_sec - t_ddMC_start.tv_sec + 1e-6*(t_ddMC_end.tv_usec - t_ddMC_start.tv_usec));
[4199b3]1210#endif 
[e2e3ad]1211                /* We loop through each row of P normalize it by making all
1212                * entries integer ones and add the resulting vector to the
1213                * int matrix facet::codim2Facets */
[2e06300]1214                for (int jj=1;jj<=/*ddakt*/P->rowsize;jj++)
[4199b3]1215                {                                       
[1d9469]1216                        fAct->numCodim2Facets++;
[4199b3]1217                        if(fAct->numCodim2Facets==1)                                   
1218                        {                                               
1219                                fAct->codim2Ptr = new facet(2);                                         
[1d9469]1220                                codim2Act = fAct->codim2Ptr;
1221                        }
1222                        else
1223                        {
1224                                codim2Act->next = new facet(2);
1225                                codim2Act = codim2Act->next;
1226                        }
[590f813]1227                        int64vec *n = new int64vec(this->numVars);
[42660f]1228#ifdef gfanp
1229                        timeval t_mI_start, t_mI_end;
1230                        gettimeofday(&t_mI_start,0);
1231#endif
[1d9469]1232                        makeInt(P,jj,*n);
[2e06300]1233                        /*for(int kk=0;kk<this->numVars;kk++)
1234                        {
1235                                int foo;
1236                                foo = (int)mpq_get_d(ddakt->matrix[ii][kk+1]);
1237                                (*n)[kk]=foo;
1238                        }*/
[42660f]1239#ifdef gfanp
1240                        gettimeofday(&t_mI_end,0);
1241                        t_mI += (t_mI_end.tv_sec - t_mI_start.tv_sec + 1e-6*(t_mI_end.tv_usec - t_mI_start.tv_usec));
1242#endif
[1d9469]1243                        codim2Act->setFacetNormal(n);
[4199b3]1244                        delete n;                                       
1245                }               
[1d9469]1246                /*We check whether the facet spanned by the codim-2 facets
1247                * intersects with the positive orthant. Otherwise we define this
[4199b3]1248                * facet to be non-flippable. Works since we set the appropriate
[2e06300]1249                * linearity for ddakt above.
[1d9469]1250                */
[4199b3]1251                //TODO It might be faster to compute jus the implied equations instead of a relative interior point
1252//              int64vec *iv_intPoint = new int64vec(this->numVars);
1253//              dd_MatrixPtr shiftMatrix;
1254//              dd_MatrixPtr intPointMatrix;
1255//              shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
1256//              for(int kk=0;kk<this->numVars;kk++)
1257//              {
1258//                      dd_set_si(shiftMatrix->matrix[kk][0],1);
1259//                      dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
1260//              }
1261//              intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
1262// #ifdef gfanp
1263//              timeval t_iP_start, t_iP_end;
1264//              gettimeofday(&t_iP_start, 0);
1265// #endif               
1266//              interiorPoint(intPointMatrix,*iv_intPoint);
1267// //           dd_rowset impl_linste,lbasis;
1268// //           dd_LPSolutionPtr lps=NULL;
1269// //           dd_ErrorType err;
1270// //           dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
1271// #ifdef gfanp
1272//              gettimeofday(&t_iP_end, 0);
1273//              t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
1274// #endif
1275//              for(int ll=0;ll<this->numVars;ll++)
1276//              {
1277//                      if( (*iv_intPoint)[ll] < 0 )
1278//                      {
1279//                              fAct->isFlippable=FALSE;
1280//                              break;
1281//                      }
1282//              }
[a7efc97]1283                /*End of check*/
1284                /*This test should be way less time consuming*/
[42660f]1285#ifdef gfanp
1286                timeval t_iP_start, t_iP_end;
1287                gettimeofday(&t_iP_start, 0);
1288#endif
[a7efc97]1289                bool containsStrictlyPosRay=TRUE;
1290                for(int ii=0;ii<ddakt->rowsize;ii++)
[1d9469]1291                {
[a7efc97]1292                        containsStrictlyPosRay=TRUE;
1293                        for(int jj=1;jj<this->numVars;jj++)
[1d9469]1294                        {
[a7efc97]1295                                if(ddakt->matrix[ii][jj]<=0)
1296                                {
1297                                        containsStrictlyPosRay=FALSE;
1298                                        break;
1299                                }
[1d9469]1300                        }
[a7efc97]1301                        if(containsStrictlyPosRay==TRUE)
1302                                break;
[f26ea0]1303                }
[a7efc97]1304                if(containsStrictlyPosRay==FALSE)
[4199b3]1305                //TODO Not sufficient. Intersect with pos orthant for pos int
[a7efc97]1306                        fAct->isFlippable=FALSE;
1307#ifdef gfanp
1308                gettimeofday(&t_iP_end, 0);
1309                t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
1310#endif
1311                /**/
[4199b3]1312                fAct = fAct->next;     
[1d9469]1313                dd_FreeMatrix(ddakt);
[6d0c34]1314                dd_FreeMatrix(P);
[4199b3]1315        }//for 
[e2e3ad]1316        dd_FreeMatrix(ddineq);
[42660f]1317#ifdef gfanp
1318        gettimeofday(&end, 0);
1319        time_getCodim2Normals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1320#endif
[1d9469]1321}
[4199b3]1322               
[72e467]1323/** Really extremal rays this time ;)
1324* Extremal rays are unique modulo the homogeneity space.
1325* Therefore we dd_MatrixAppend gc->ddFacets and gcone::dd_LinealitySpace
1326* into ddineq. Next we compute the extremal rays of the so given subspace.
1327* Figuring out whether a ray belongs to a given facet(normal) is done by
1328* checking whether the inner product of the ray with the normal is zero.
[590f813]1329* We use ivAdd here which returns a new int64vec. Therefore we need to avoid
[4199b3]1330* a memory leak which would be cause by the line
[2f72c55]1331* iv=ivAdd(iv,b)
[4199b3]1332* So we keep pointer tmp to iv and delete(tmp), so there should not occur a
[2f72c55]1333* memleak
[73809b]1334* TODO normalization
[72e467]1335*/
1336void gcone::getExtremalRays(const gcone &gc)
1337{
[4c4a80]1338#ifdef gfanp
1339        timeval start, end;
1340        gettimeofday(&start, 0);
[73809b]1341        timeval poly_start, poly_end;
1342        gettimeofday(&poly_start,0);
[4c4a80]1343#endif
[4199b3]1344        //Add lineality space - dd_LinealitySpace
[72e467]1345        dd_MatrixPtr ddineq;
[4199b3]1346        dd_ErrorType err;       
[10dd14]1347        ddineq = (dd_LinealitySpace->rowsize>0) ? dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace) : dd_CopyMatrix(gc.ddFacets);
[be7ab3f]1348        /* In case the input is non-homogeneous we add constrains for the positive orthant.
[4199b3]1349        * This is justified by the fact that for non-homog ideals we only consider the
[be7ab3f]1350        * restricted fan. This way we can be sure to find strictly positive interior points.
1351        * This in turn makes life easy when checking for flippability!
1352        * Drawback: Makes the LP larger so probably slows down computations a wee bit.
1353        */
1354        dd_MatrixPtr ddPosRestr;
1355        if(hasHomInput==FALSE)
1356        {
1357                dd_MatrixPtr tmp;
1358                ddPosRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1359                for(int ii=0;ii<this->numVars;ii++)
1360                        dd_set_si(ddPosRestr->matrix[ii][ii+1],1);
1361                dd_MatrixAppendTo(&ddineq,ddPosRestr);
1362                assert(ddineq);
1363                dd_FreeMatrix(ddPosRestr);
[4199b3]1364        }       
[72e467]1365        dd_PolyhedraPtr ddPolyh;
1366        ddPolyh = dd_DDMatrix2Poly(ddineq, &err);
1367        dd_MatrixPtr P;
[e58dbb]1368        P=dd_CopyGenerators(ddPolyh);//Here we actually compute the rays!
[be7ab3f]1369        dd_FreePolyhedra(ddPolyh);
[73809b]1370        dd_FreeMatrix(ddineq);
1371#ifdef gfanp
1372        gettimeofday(&poly_end,0);
1373        t_ddPolyh += (poly_end.tv_sec - poly_start.tv_sec + 1e-6*(poly_end.tv_usec - poly_start.tv_usec));
1374#endif
[be7ab3f]1375        /* Compute interior point on the fly*/
[590f813]1376        int64vec *ivIntPointOfCone = new int64vec(this->numVars);
[dd6b1d]1377        int64vec *foo = new int64vec(this->numVars);
[590f813]1378        for(int ii=0;ii<P->rowsize;ii++)
1379        {
1380                int64vec *tmp = ivIntPointOfCone;
1381                makeInt(P,ii+1,*foo);
1382                ivIntPointOfCone = iv64Add(ivIntPointOfCone,foo);
1383                delete tmp;
1384        }
[dd6b1d]1385        delete foo;
[590f813]1386        int64 ggT=(*ivIntPointOfCone)[0];
[be7ab3f]1387        for (int ii=0;ii<(this->numVars);ii++)
1388        {
[4199b3]1389                if( (*ivIntPointOfCone)[ii]>INT_MAX ) 
[590f813]1390                        WarnS("Interior point exceeds INT_MAX!\n");
[4199b3]1391                //Compute intgcd
[f91fddc]1392                ggT=int64gcd(ggT,(*ivIntPointOfCone)[ii]);
[73809b]1393        }
[4199b3]1394       
1395        //Divide out a common gcd > 1
[73809b]1396        if(ggT>1)
1397        {
1398                for(int ii=0;ii<this->numVars;ii++)
[590f813]1399                {
[73809b]1400                        (*ivIntPointOfCone)[ii] /= ggT;
[590f813]1401                        if( (*ivIntPointOfCone)[ii]>INT_MAX ) WarnS("Interior point still exceeds INT_MAX after GCD!\n");
1402                }
[be7ab3f]1403        }
[e58dbb]1404
[be7ab3f]1405        /*For homogeneous input (like Det3,3,5) the int points may be negative. So add a suitable multiple of (1,_,1)*/
1406        if(hasHomInput==TRUE && iv64isStrictlyPositive(ivIntPointOfCone)==FALSE)
1407        {
[590f813]1408                int64vec *ivOne = new int64vec(this->numVars);
[ec6c52]1409                int maxNegEntry=0;
[be7ab3f]1410                for(int ii=0;ii<this->numVars;ii++)
1411                {
[4199b3]1412//                      (*ivOne)[ii]=1;
[ec6c52]1413                        if ((*ivIntPointOfCone)[ii]<maxNegEntry) maxNegEntry=(*ivIntPointOfCone)[ii];
[be7ab3f]1414                }
[ec6c52]1415                maxNegEntry *= -1;
[4199b3]1416                maxNegEntry++;//To be on the safe side
[ec6c52]1417                for(int ii=0;ii<this->numVars;ii++)
1418                        (*ivOne)[ii]=maxNegEntry;
[590f813]1419                int64vec *tmp=ivIntPointOfCone;
1420                ivIntPointOfCone=iv64Add(ivIntPointOfCone,ivOne);
[ec6c52]1421                delete(tmp);
[4199b3]1422//              while( !iv64isStrictlyPositive(ivIntPointOfCone) )
1423//              {
1424//                      int64vec *tmp = ivIntPointOfCone;
1425//                      for(int jj=0;jj<this->numVars;jj++)
1426//                              (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
1427//                      ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
1428//                      delete tmp;                             
1429//              }
[be7ab3f]1430                delete ivOne;
[590f813]1431                int64 ggT=(*ivIntPointOfCone)[0];
[73809b]1432                for(int ii=0;ii<this->numVars;ii++)
[f91fddc]1433                        ggT=int64gcd( ggT, (*ivIntPointOfCone)[ii]);
[73809b]1434                if(ggT>1)
1435                {
1436                        for(int jj=0;jj<this->numVars;jj++)
1437                                (*ivIntPointOfCone)[jj] /= ggT;
1438                }
[be7ab3f]1439        }
[4199b3]1440//      assert(iv64isStrictlyPositive(ivIntPointOfCone));
1441       
[be7ab3f]1442        this->setIntPoint(ivIntPointOfCone);
1443        delete(ivIntPointOfCone);
1444        /* end of interior point computation*/
[4199b3]1445       
1446        //Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
[72e467]1447        int rows=P->rowsize;
1448        facet *fAct=gc.facetPtr;
[4199b3]1449        //Construct an array to hold the extremal rays of the cone
1450        this->gcRays = (int64vec**)omAlloc0(sizeof(int64vec*)*P->rowsize);     
[68eb0c]1451        for(int ii=0;ii<P->rowsize;ii++)
1452        {
[590f813]1453                int64vec *rowvec = new int64vec(this->numVars);
[4199b3]1454                makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve
[590f813]1455                this->gcRays[ii] = iv64Copy(rowvec);
[68eb0c]1456                delete rowvec;
[4199b3]1457        }       
[68eb0c]1458        this->numRays=P->rowsize;
[4199b3]1459        //Check which rays belong to which facet
[72e467]1460        while(fAct!=NULL)
1461        {
[4199b3]1462                const int64vec *fNormal;// = new int64vec(this->numVars);
1463                fNormal = fAct->getRef2FacetNormal();//->getFacetNormal();
[590f813]1464                int64vec *ivIntPointOfFacet = new int64vec(this->numVars);
[72e467]1465                for(int ii=0;ii<rows;ii++)
[4199b3]1466                {                       
[68eb0c]1467                        if(dotProduct(*fNormal,this->gcRays[ii])==0)
[72e467]1468                        {
[4199b3]1469                                int64vec *tmp = ivIntPointOfFacet;//Prevent memleak
[72e467]1470                                fAct->numCodim2Facets++;
1471                                facet *codim2Act;
[4199b3]1472                                if(fAct->numCodim2Facets==1)                                   
1473                                {                                               
1474                                        fAct->codim2Ptr = new facet(2);                                         
[72e467]1475                                        codim2Act = fAct->codim2Ptr;
1476                                }
1477                                else
1478                                {
1479                                        codim2Act->next = new facet(2);
1480                                        codim2Act = codim2Act->next;
1481                                }
[4199b3]1482                                //codim2Act->setFacetNormal(rowvec);
1483                                //Rather just let codim2Act point to the corresponding int64vec of gcRays
[68eb0c]1484                                codim2Act->fNormal=this->gcRays[ii];
[4199b3]1485                                fAct->numRays++;                                 
1486                                //Memleak avoided via tmp
[590f813]1487                                ivIntPointOfFacet=iv64Add(ivIntPointOfFacet,this->gcRays[ii]);
[4199b3]1488                                //Now tmp still points to the OLD address of ivIntPointOfFacet
[2f72c55]1489                                delete(tmp);
[4199b3]1490                                       
[72e467]1491                        }
[4199b3]1492                }//For non-homog input ivIntPointOfFacet should already be >0 here
1493//              if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
1494                //if we have no strictly pos ray but the input is homogeneous
1495                //then add a suitable multiple of (1,...,1)
[be7ab3f]1496                if( !iv64isStrictlyPositive(ivIntPointOfFacet) && hasHomInput==TRUE)
1497                {
[590f813]1498                        int64vec *ivOne = new int64vec(this->numVars);
[be7ab3f]1499                        for(int ii=0;ii<this->numVars;ii++)
1500                                (*ivOne)[ii]=1;
1501                        while( !iv64isStrictlyPositive(ivIntPointOfFacet) )
1502                        {
[590f813]1503                                int64vec *tmp = ivIntPointOfFacet;
[be7ab3f]1504                                for(int jj=0;jj<this->numVars;jj++)
1505                                {
[4199b3]1506                                        (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
[be7ab3f]1507                                }
[590f813]1508                                ivIntPointOfFacet = iv64Add(ivIntPointOfFacet/*diff*/,ivOne);
[4199b3]1509                                delete tmp;                             
[be7ab3f]1510                        }
1511                        delete ivOne;
1512                }
[590f813]1513                int64 ggT=(*ivIntPointOfFacet)[0];
[73809b]1514                for(int ii=0;ii<this->numVars;ii++)
[f91fddc]1515                        ggT=int64gcd(ggT,(*ivIntPointOfFacet)[ii]);
[73809b]1516                if(ggT>1)
1517                {
1518                        for(int ii=0;ii<this->numVars;ii++)
1519                                (*ivIntPointOfFacet)[ii] /= ggT;
[4199b3]1520                }                       
[73809b]1521                fAct->setInteriorPoint(ivIntPointOfFacet);
[4199b3]1522               
[2f72c55]1523                delete(ivIntPointOfFacet);
[4199b3]1524                //Now (if we have at least 3 variables) do a bubblesort on the rays
[5ece85]1525                /*if(this->numVars>2)
1526                {
1527                        facet *A[fAct->numRays-1];
1528                        facet *f2Act=fAct->codim2Ptr;
1529                        for(unsigned ii=0;ii<fAct->numRays;ii++)
1530                        {
1531                                A[ii]=f2Act;
1532                                f2Act=f2Act->next;
1533                        }
1534                        bool exchanged=FALSE;
1535                        unsigned n=fAct->numRays-1;
1536                        do
1537                        {
[4199b3]1538                                exchanged=FALSE;//n=fAct->numRays-1;                           
[5ece85]1539                                for(unsigned ii=0;ii<=n-1;ii++)
[4199b3]1540                                {                                       
[5ece85]1541                                        if((A[ii]->fNormal)->compare((A[ii+1]->fNormal))==1)
1542                                        {
[4199b3]1543                                                //Swap rays
[5ece85]1544                                                cout << "Swapping ";
1545                                                A[ii]->fNormal->show(1,0); cout << " with "; A[ii+1]->fNormal->show(1,0); cout << endl;
1546                                                A[ii]->next=A[ii+1]->next;
1547                                                if(ii>0)
1548                                                        A[ii-1]->next=A[ii+1];
1549                                                A[ii+1]->next=A[ii];
1550                                                if(ii==0)
1551                                                        fAct->codim2Ptr=A[ii+1];
[4199b3]1552                                                //end swap
1553                                                facet *tmp=A[ii];//swap in list
[5ece85]1554                                                A[ii+1]=A[ii];
1555                                                A[ii]=tmp;
[4199b3]1556//                                              tmp=NULL;                       
1557                                        }                                       
[5ece85]1558                                }
[4199b3]1559                                n--;                   
[5ece85]1560                        }while(exchanged==TRUE && n>=0);
[1f637e]1561                }*///if (currRing->N)>2
[4199b3]1562//              delete fNormal;         
[72e467]1563                fAct = fAct->next;
[4199b3]1564        }//end of facet checking
[be7ab3f]1565        dd_FreeMatrix(P);
[4199b3]1566        //Now all extremal rays should be set w.r.t their respective fNormal
1567        //TODO Not sufficient -> vol2 II/125&127
1568        //NOTE Sufficient according to cddlibs doc. These ARE rays
[871f01]1569        //What the hell... let's just take interior points
[85583b]1570        if(gcone::hasHomInput==FALSE)
1571        {
[5e4025]1572                fAct=gc.facetPtr;
1573                while(fAct!=NULL)
[72e467]1574                {
[4199b3]1575//                      bool containsStrictlyPosRay=FALSE;
1576//                      facet *codim2Act;
1577//                      codim2Act = fAct->codim2Ptr;
1578//                      while(codim2Act!=NULL)
1579//                      {
1580//                              int64vec *rayvec;
1581//                              rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
1582//                              //int negCtr=0;
1583//                              if(iv64isStrictlyPositive(rayvec))
1584//                              {
1585//                                      containsStrictlyPosRay=TRUE;
1586//                                      delete(rayvec);
1587//                                      break;
1588//                              }                               
1589//                              delete(rayvec);
1590//                              codim2Act = codim2Act->next;
1591//                      }
1592//                      if(containsStrictlyPosRay==FALSE)
1593//                              fAct->isFlippable=FALSE;
[be7ab3f]1594                        if(!iv64isStrictlyPositive(fAct->interiorPoint))
[5e4025]1595                                fAct->isFlippable=FALSE;
1596                        fAct = fAct->next;
[72e467]1597                }
[4199b3]1598        }//hasHomInput?
[4c4a80]1599#ifdef gfanp
1600        gettimeofday(&end, 0);
1601        t_getExtremalRays += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
[4199b3]1602#endif 
[305057]1603}
1604
1605/** Order the spanning rays in a lex way hopefully using qsort()*/
1606void gcone::orderRays()
1607{
[4199b3]1608//   qsort(gcRays,sizeof(int64vec),int64vec::compare);
[305057]1609}
[4199b3]1610       
[590f813]1611inline bool gcone::iv64isStrictlyPositive(const int64vec * iv64)
[85583b]1612{
1613        bool res=TRUE;
1614        for(int ii=0;ii<iv64->length();ii++)
1615        {
1616                if((*iv64)[ii]<=0)
1617                {
1618                        res=FALSE;
1619                        break;
[4199b3]1620                }               
[85583b]1621        }
1622        return res;
1623}
[305057]1624
[4199b3]1625/** \brief Compute the Groebner Basis on the other side of a shared facet
[1d9469]1626 *
1627 * Implements algorithm 4.3.2 from Anders' thesis.
1628 * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
1629 * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
1630 * is parallel to \f$ leadexp(g) \f$
1631 * Parallelity is checked using basic linear algebra. See gcone::isParallel.
1632 * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
[4199b3]1633 * computing an interior point of the facet and taking all terms having the same weight with respect
[1d9469]1634 * to this interior point.
1635 *\param ideal, facet
1636 * Input: a marked,reduced Groebner basis and a facet
1637 */
[4199b3]1638inline void gcone::flip(ideal gb, facet *f)             //Compute "the other side"
1639{       
[42660f]1640#ifdef gfanp
1641        timeval start, end;
1642        gettimeofday(&start, 0);
[4199b3]1643#endif         
1644        int64vec *fNormal;// = new int64vec(this->numVars);     //facet normal, check for parallelity                   
1645        fNormal = f->getFacetNormal();  //read this->fNormal;
[96e408]1646#ifndef NDEBUG
[4199b3]1647//      std::cout << "running gcone::flip" << std::endl;
[11a7dc]1648        printf("flipping UCN %i over facet",this->getUCN());
[be7ab3f]1649        fNormal->show(1,0);
[4199b3]1650        printf(") with UCN %i\n",f->getUCN() ); 
[be7ab3f]1651#endif
[1d37196]1652        if(this->getUCN() != f->getUCN())
1653        {
1654                WerrorS("Uh oh... Trying to flip over facet with incompatible UCN");
1655                exit(-1);
1656        }
[1d9469]1657        /*1st step: Compute the initial ideal*/
[4199b3]1658        /*poly initialFormElement[IDELEMS(gb)];*/       //array of #polys in GB to store initial form
[1d9469]1659        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
[4199b3]1660       
[7c0ec6]1661        computeInv(gb,initialForm,*fNormal);
[1d37196]1662
[96e408]1663#ifndef NDEBUG
[bb503c7]1664/*      cout << "Initial ideal is: " << endl;
[1d9469]1665        idShow(initialForm);
[4199b3]1666        //f->printFlipGB();*/
1667//      cout << "===" << endl;
1668#endif                 
[1d9469]1669        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
[bb503c7]1670        /*Substep 2.1
[4199b3]1671        compute $G_{-\alpha}(in_v(I))
[1d9469]1672        see journal p. 66
1673        NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
1674        srcRing already has a weighted ordering
[bb503c7]1675        */
[1d9469]1676        ring srcRing=currRing;
1677        ring tmpRing;
[4199b3]1678                       
[1d9469]1679        if( (srcRing->order[0]!=ringorder_a))
1680        {
[4199b3]1681                int64vec *iv;// = new int64vec(this->numVars);
1682                iv = ivNeg(fNormal);//ivNeg uses iv64Copy -> new
1683//              tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
[dcf8b4b]1684                tmpRing=rCopyAndAddWeight(srcRing,iv);
1685                delete iv;
[1d9469]1686        }
1687        else
1688        {
1689                tmpRing=rCopy0(srcRing);
1690                int length=fNormal->length();
1691                int *A=(int *)omAlloc0(length*sizeof(int));
1692                for(int jj=0;jj<length;jj++)
1693                {
1694                        A[jj]=-(*fNormal)[jj];
1695                }
1696                omFree(tmpRing->wvhdl[0]);
1697                tmpRing->wvhdl[0]=(int*)A;
1698                tmpRing->block1[0]=length;
[a23a297]1699                rComplete(tmpRing);
[4199b3]1700                //omFree(A);
[1d9469]1701        }
[4199b3]1702        delete fNormal; 
1703        rChangeCurrRing(tmpRing);       
1704                       
1705        ideal ina;                     
[7c0ec6]1706        ina=idrCopyR(initialForm,srcRing);
[1d37196]1707        idDelete(&initialForm);
[1d9469]1708        ideal H;
[4199b3]1709//      H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
[72e467]1710#ifdef gfanp
1711        timeval t_kStd_start, t_kStd_end;
1712        gettimeofday(&t_kStd_start,0);
1713#endif
[2f72c55]1714        if(gcone::hasHomInput==TRUE)
1715                H=kStd(ina,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
1716        else
[4199b3]1717                H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
[72e467]1718#ifdef gfanp
1719        gettimeofday(&t_kStd_end, 0);
1720        t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
1721#endif
[1d9469]1722        idSkipZeroes(H);
[1d37196]1723        idDelete(&ina);
1724
[bb503c7]1725        /*Substep 2.2
[1d9469]1726        do the lifting and mark according to H
[bb503c7]1727        */
[1d9469]1728        rChangeCurrRing(srcRing);
1729        ideal srcRing_H;
[4199b3]1730        ideal srcRing_HH;                       
[1d9469]1731        srcRing_H=idrCopyR(H,tmpRing);
[871f01]1732        //H is needed further below, so don't idDelete here
[1d37196]1733        srcRing_HH=ffG(srcRing_H,this->gcBasis);
[dcf8b4b]1734        idDelete(&srcRing_H);
[4199b3]1735               
[bb503c7]1736        /*Substep 2.2.1
1737         * Mark according to G_-\alpha
1738         * Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
1739         * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
1740         * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
[4199b3]1741         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
[bb503c7]1742         * compute the difference accordingly
1743        */
[2e06300]1744#ifdef gfanp
1745        timeval t_markings_start, t_markings_end;
1746        gettimeofday(&t_markings_start, 0);
[4199b3]1747#endif         
[1d9469]1748        bool markingsAreCorrect=FALSE;
1749        dd_MatrixPtr intPointMatrix;
1750        int iPMatrixRows=0;
[4199b3]1751        dd_rowrange aktrow=0;                   
[1d9469]1752        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1753        {
[871f01]1754                poly aktpoly=(poly)srcRing_HH->m[ii];//This is a pointer, so don't pDelete
[4199b3]1755                iPMatrixRows = iPMatrixRows+pLength(aktpoly);           
[1d9469]1756        }
[bb503c7]1757        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
1758         * construction of the differences
1759         */
[4199b3]1760        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); 
1761        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
1762       
[1d9469]1763        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1764        {
[4199b3]1765                markingsAreCorrect=FALSE;       //crucial to initialise here
[871f01]1766                poly aktpoly=srcRing_HH->m[ii]; //Only a pointer, so don't pDelete
[1d9469]1767                /*Comparison of leading monomials is done via exponent vectors*/
1768                for (int jj=0;jj<IDELEMS(H);jj++)
1769                {
1770                        int *src_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
1771                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
1772                        pGetExpV(aktpoly,src_ExpV);
[4199b3]1773                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
[f91fddc]1774                        poly p=pCopy(H->m[ii]);
1775                        pGetExpV(p/*pCopy(H->m[ii])*/,dst_ExpV);
1776                        pDelete(&p);
[1d9469]1777                        rChangeCurrRing(srcRing);
1778                        bool expVAreEqual=TRUE;
1779                        for (int kk=1;kk<=this->numVars;kk++)
[459b80]1780                        {
[96e408]1781#ifndef NDEBUG
[4199b3]1782//                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
[348034]1783#endif
[1d9469]1784                                if (src_ExpV[kk]!=dst_ExpV[kk])
1785                                {
1786                                        expVAreEqual=FALSE;
1787                                }
[4199b3]1788                        }                                                               
[1d9469]1789                        if (expVAreEqual==TRUE)
1790                        {
[4199b3]1791                                markingsAreCorrect=TRUE; //everything is fine
[96e408]1792#ifndef NDEBUG
[4199b3]1793//                              cout << "correct markings" << endl;
[348034]1794#endif
[4199b3]1795                        }//if (pHead(aktpoly)==pHead(H->m[jj])
[7c0ec6]1796                        omFree(src_ExpV);
1797                        omFree(dst_ExpV);
[4199b3]1798                }//for (int jj=0;jj<IDELEMS(H);jj++)
1799               
[1d9469]1800                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
1801                if (markingsAreCorrect==TRUE)
1802                {
1803                        pGetExpV(aktpoly,leadExpV);
1804                }
1805                else
1806                {
1807                        rChangeCurrRing(tmpRing);
[4199b3]1808                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
[1d9469]1809                        rChangeCurrRing(srcRing);
1810                }
[4199b3]1811                /*compute differences of the expvects*/                         
[1d9469]1812                while (pNext(aktpoly)!=NULL)
1813                {
[42660f]1814                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
[4199b3]1815                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
[1d9469]1816                        is not omitted when computing the differences*/
1817                        if(markingsAreCorrect==TRUE)
[b1bed35]1818                        {
[1d9469]1819                                aktpoly=pNext(aktpoly);
1820                                pGetExpV(aktpoly,v);
[b1bed35]1821                        }
[1d9469]1822                        else
[4e4d75]1823                        {
[1d9469]1824                                pGetExpV(pHead(aktpoly),v);
1825                                markingsAreCorrect=TRUE;
[4e4d75]1826                        }
[a7efc97]1827                        int ctr=0;
[1d9469]1828                        for (int jj=0;jj<this->numVars;jj++)
[4199b3]1829                        {                               
1830                                /*Store into ddMatrix*/                         
[a7efc97]1831                                if(leadExpV[jj+1]-v[jj+1])
1832                                        ctr++;
[1d9469]1833                                dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
1834                        }
[a7efc97]1835                        /*It ought to be more sensible to avoid 0-rows in the first place*/
[4199b3]1836//                      if(ctr==this->numVars)//We have a 0-row
1837//                              dd_MatrixRowRemove(&intPointMatrix,aktrow);
1838//                      else
[a7efc97]1839                                aktrow +=1;
[42660f]1840                        omFree(v);
[1d9469]1841                }
[7c0ec6]1842                omFree(leadExpV);
[4199b3]1843        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
[2e06300]1844#ifdef gfanp
1845        gettimeofday(&t_markings_end, 0);
1846        t_markings += (t_markings_end.tv_sec - t_markings_start.tv_sec + 1e-6*(t_markings_end.tv_usec - t_markings_start.tv_usec));
1847#endif
[1d37196]1848        /*Now it is safe to idDelete(H)*/
1849        idDelete(&H);
[a7efc97]1850        /*Preprocessing goes here since otherwise we would delete the constraint
1851        * for the standard simplex.
1852        */
1853        preprocessInequalities(intPointMatrix);
[2e06300]1854        /*Now we add the constraint for the standard simplex*/
[4199b3]1855//      dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
1856//      for (int jj=1;jj<=this->numVars;jj++)
1857//      {
1858//              dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
1859//      }
[871f01]1860        //Let's make sure we compute interior points from the positive orthant
[4199b3]1861//      dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1862//     
1863//      int jj=1;
1864//      for (int ii=0;ii<this->numVars;ii++)
1865//      {
1866//              dd_set_si(posRestr->matrix[ii][jj],1);
1867//              jj++;                                                   
1868//      }
[a7efc97]1869        /*We create a matrix containing the standard simplex
1870        * and constraints to assure a strictly positive point
[72e467]1871        * is computed */
[a7efc97]1872        dd_MatrixPtr posRestr = dd_CreateMatrix(this->numVars+1, this->numVars+1);
1873        for(int ii=0;ii<posRestr->rowsize;ii++)
[1d9469]1874        {
[a7efc97]1875                if(ii==0)
1876                {
1877                        dd_set_si(posRestr->matrix[ii][0],-1);
[4199b3]1878                        for(int jj=1;jj<=this->numVars;jj++)                   
1879                                dd_set_si(posRestr->matrix[ii][jj],1);                 
[a7efc97]1880                }
1881                else
1882                {
[72e467]1883                        /** Set all variables to \geq 1/10. YMMV but this choice is pretty equal*/
1884                        dd_set_si2(posRestr->matrix[ii][0],-1,2);
[a7efc97]1885                        dd_set_si(posRestr->matrix[ii][ii],1);
1886                }
[1d9469]1887        }
1888        dd_MatrixAppendTo(&intPointMatrix,posRestr);
1889        dd_FreeMatrix(posRestr);
[72e467]1890
[590f813]1891        int64vec *iv_weight = new int64vec(this->numVars);
[2e06300]1892#ifdef gfanp
1893        timeval t_dd_start, t_dd_end;
1894        gettimeofday(&t_dd_start, 0);
1895#endif
1896        dd_ErrorType err;
1897        dd_rowset implLin, redrows;
1898        dd_rowindex newpos;
[72e467]1899
[4199b3]1900        //NOTE Here we should remove interiorPoint and instead
1901        // create and ordering like (a(omega),a(fNormal),dp)
1902//      if(this->ivIntPt==NULL)
1903                interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
1904//      else
1905//              iv_weight=this->getIntPoint();
[1d9469]1906        dd_FreeMatrix(intPointMatrix);
[2e06300]1907        /*Crude attempt for interior point */
1908        /*dd_PolyhedraPtr ddpolyh;
1909        dd_ErrorType err;
1910        dd_rowset impl_linset,redset;
1911        dd_rowindex newpos;
1912        dd_MatrixCanonicalize(&intPointMatrix, &impl_linset, &redset, &newpos, &err);
1913        ddpolyh=dd_DDMatrix2Poly(intPointMatrix, &err);
1914        dd_MatrixPtr P;
1915        P=dd_CopyGenerators(ddpolyh);
1916        dd_FreePolyhedra(ddpolyh);
1917        for(int ii=0;ii<P->rowsize;ii++)
1918        {
[590f813]1919                int64vec *iv_row=new int64vec(this->numVars);
[2e06300]1920                makeInt(P,ii+1,*iv_row);
1921                iv_weight =ivAdd(iv_weight, iv_row);
1922                delete iv_row;
1923        }
1924        dd_FreeMatrix(P);
1925        dd_FreeMatrix(intPointMatrix);*/
1926#ifdef gfanp
1927        gettimeofday(&t_dd_end, 0);
1928        t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec));
[4199b3]1929#endif 
1930       
[963eeb]1931        /*Step 3
[4199b3]1932        * turn the minimal basis into a reduced one */                 
1933        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
1934        // Thus:
1935        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
[1d9469]1936        ring dstRing=rCopy0(tmpRing);
1937        int length=iv_weight->length();
1938        int *A=(int *)omAlloc0(length*sizeof(int));
1939        for(int jj=0;jj<length;jj++)
1940        {
1941                A[jj]=(*iv_weight)[jj];
1942        }
1943        dstRing->wvhdl[0]=(int*)A;
[1d37196]1944        rComplete(dstRing);
[1d9469]1945        rChangeCurrRing(dstRing);
[1d37196]1946        rDelete(tmpRing);
[1d9469]1947        delete iv_weight;
[e2e3ad]1948
[4199b3]1949        ideal dstRing_I;                       
[1d9469]1950        dstRing_I=idrCopyR(srcRing_HH,srcRing);
[4199b3]1951        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more
1952        //dstRing_I=idrCopyR(inputIdeal,srcRing);
[1d9469]1953        BITSET save=test;
1954        test|=Sy_bit(OPT_REDSB);
1955        test|=Sy_bit(OPT_REDTAIL);
[96e408]1956#ifndef NDEBUG
[4199b3]1957//      test|=Sy_bit(6);        //OPT_DEBUG
[f72928]1958#endif
[1d9469]1959        ideal tmpI;
[4199b3]1960        //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
1961        //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
[1d9469]1962        tmpI = dstRing_I;
[72e467]1963#ifdef gfanp
1964        gettimeofday(&t_kStd_start,0);
1965#endif
[2f72c55]1966        if(gcone::hasHomInput==TRUE)
1967                dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
1968        else
1969                dstRing_I=kStd(tmpI,NULL,isNotHomog,NULL);
[72e467]1970#ifdef gfanp
1971        gettimeofday(&t_kStd_end, 0);
1972        t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
1973#endif
[1d37196]1974        idDelete(&tmpI);
[4199b3]1975        idNorm(dstRing_I);                     
1976//      kInterRed(dstRing_I);
[1d9469]1977        idSkipZeroes(dstRing_I);
1978        test=save;
1979        /*End of step 3 - reduction*/
[4199b3]1980                       
1981        f->setFlipGB(dstRing_I);//store the flipped GB
1982//      idDelete(&dstRing_I);
1983        f->flipRing=rCopy(dstRing);     //store the ring on the other side
[96e408]1984#ifndef NDEBUG
[11a7dc]1985        printf("Flipped GB is UCN %i:\n",counter+1);
1986        idDebugPrint(dstRing_I);
1987        printf("\n");
[4199b3]1988#endif 
1989        idDelete(&dstRing_I);   
1990        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
[1d37196]1991        rDelete(dstRing);
[42660f]1992#ifdef gfanp
1993        gettimeofday(&end, 0);
1994        time_flip += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1995#endif
[4199b3]1996}//void flip(ideal gb, facet *f)
[7c0ec6]1997
[be7ab3f]1998/** \brief A slightly different approach to flipping
1999* Here we use the fact that in_v(in_u(I))=in_(u+eps*v)(I). Therefore, we do no longer
2000* need to compute an interior point and run BBA on the minimal basis but we can rather
2001* use the ordering (a(omega),a(fNormal),dp)
2002* The second parameter facet *f must not be const since we need to store f->flipGB
2003* Problem: Assume we start in a cone with ordering (dp,C). Then \f$ in_\omega(I) \f$
2004* will be from a ring with (a(),dp,C) and our resulting cone from (a(),a(),dp,C). Hence a way
2005* must be found to circumvent the sequence of a()'s growing to a ridiculous size.
[4199b3]2006* Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we
2007* do have an interior point of the cone by adding the extremal rays. So we replace
2008* the latter cone by a cone with (a(sum_of_rays),dp,C).
[be7ab3f]2009* Con: It's incredibly ugly
2010* Pro: No messing around with readConeFromFile()
2011* Is there a way to construct a vector from \f$ \omega \f$ and the facet normal?
2012*/
[11a7dc]2013inline void gcone::flip2(const ideal &gb, facet *f)
[be7ab3f]2014{
2015#ifdef gfanp
2016        timeval start, end;
2017        gettimeofday(&start, 0);
2018#endif
[590f813]2019        const int64vec *fNormal;
[4199b3]2020        fNormal = f->getRef2FacetNormal();/*->getFacetNormal();*/       //read this->fNormal;
[96e408]2021#ifndef NDEBUG
[d503ac]2022//      printf("flipping UCN %i over facet(",this->getUCN());
2023//      fNormal->show(1,0);
2024//      printf(") with UCN %i\n",f->getUCN()); 
[be7ab3f]2025#endif
2026        if(this->getUCN() != f->getUCN())
[11a7dc]2027        {       printf("%i vs %i\n",this->getUCN(), f->getUCN() );
[be7ab3f]2028                WerrorS("Uh oh... Trying to flip over facet with incompatible UCN");
2029                exit(-1);
2030        }
2031        /*1st step: Compute the initial ideal*/
[4199b3]2032        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);     
[be7ab3f]2033        computeInv( gb, initialForm, *fNormal );
2034        ring srcRing=currRing;
2035        ring tmpRing;
[4199b3]2036       
[590f813]2037        const int64vec *intPointOfFacet;
[be7ab3f]2038        intPointOfFacet=f->getInteriorPoint();
[4199b3]2039        //Now we need two blocks of ringorder_a!
2040        //May assume the same situation as in flip() here                       
[590f813]2041        if( (srcRing->order[0]!=ringorder_a/*64*/) && (srcRing->order[1]!=ringorder_a/*64*/) )
[be7ab3f]2042        {
[4199b3]2043                int64vec *iv = new int64vec(this->numVars);//init with 1s, since we do not need a 2nd block here but later
2044//              int64vec *iv_foo = new int64vec(this->numVars,1);//placeholder
2045                int64vec *ivw = ivNeg(const_cast<int64vec*>(fNormal));         
[be7ab3f]2046                tmpRing=rCopyAndAddWeight2(srcRing,ivw/*intPointOfFacet*/,iv);
2047                delete iv;delete ivw;
[4199b3]2048//              delete iv_foo;
[be7ab3f]2049        }
[4199b3]2050        else 
[be7ab3f]2051        {
[590f813]2052                int64vec *iv=new int64vec(this->numVars);
2053                int64vec *ivw=ivNeg(const_cast<int64vec*>(fNormal));
[be7ab3f]2054                tmpRing=rCopyAndAddWeight2(srcRing,ivw,iv);
2055                delete iv; delete ivw;
2056                /*tmpRing=rCopy0(srcRing);
2057                int length=fNormal->length();
2058                int *A1=(int *)omAlloc0(length*sizeof(int));
2059                int *A2=(int *)omAlloc0(length*sizeof(int));
2060                for(int jj=0;jj<length;jj++)
2061                {
2062                        A1[jj] = -(*fNormal)[jj];
[4199b3]2063                        A2[jj] = 1;//-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
[be7ab3f]2064                }
2065                omFree(tmpRing->wvhdl[0]);
2066                if(tmpRing->wvhdl[1]!=NULL)
2067                        omFree(tmpRing->wvhdl[1]);
[4199b3]2068                tmpRing->wvhdl[0]=(int*)A1;             
[be7ab3f]2069                tmpRing->block1[0]=length;
[4199b3]2070                tmpRing->wvhdl[1]=(int*)A2;             
[be7ab3f]2071                tmpRing->block1[1]=length;
2072                rComplete(tmpRing);*/
2073        }
[4199b3]2074//      delete fNormal; //NOTE Do not delete when using getRef2FacetNormal();
2075        rChangeCurrRing(tmpRing);       
2076        //Now currRing should have (a(),a(),dp,C)               
2077        ideal ina;                     
[be7ab3f]2078        ina=idrCopyR(initialForm,srcRing);
2079        idDelete(&initialForm);
2080        ideal H;
2081#ifdef gfanp
2082        timeval t_kStd_start, t_kStd_end;
2083        gettimeofday(&t_kStd_start,0);
2084#endif
2085        BITSET save=test;
2086        test|=Sy_bit(OPT_REDSB);
2087        test|=Sy_bit(OPT_REDTAIL);
[4199b3]2088//      if(gcone::hasHomInput==TRUE)
[be7ab3f]2089                H=kStd(ina,NULL,testHomog/*isHomog*/,NULL/*,gcone::hilbertFunction*/);
[4199b3]2090//      else
2091//              H=kStd(ina,NULL,isNotHomog,NULL);       //This is \mathcal(G)_{>_-\alpha}(in_v(I))
[be7ab3f]2092        test=save;
2093#ifdef gfanp
2094        gettimeofday(&t_kStd_end, 0);
2095        t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
2096#endif
2097        idSkipZeroes(H);
2098        idDelete(&ina);
[4199b3]2099       
[be7ab3f]2100        rChangeCurrRing(srcRing);
2101        ideal srcRing_H;
[4199b3]2102        ideal srcRing_HH;                       
[be7ab3f]2103        srcRing_H=idrCopyR(H,tmpRing);
[871f01]2104        //H is needed further below, so don't idDelete here
[be7ab3f]2105        srcRing_HH=ffG(srcRing_H,this->gcBasis);
2106        idDelete(&srcRing_H);
[4199b3]2107        //Now BBA(srcRing_HH) with (a(),a(),dp)
[be7ab3f]2108        /* Evil modification of currRing */
2109        ring dstRing=rCopy0(tmpRing);
2110        int length=this->numVars;
2111        int *A1=(int *)omAlloc0(length*sizeof(int));
2112        int *A2=(int *)omAlloc0(length*sizeof(int));
[590f813]2113        const int64vec *ivw=f->getRef2FacetNormal();
[be7ab3f]2114        for(int jj=0;jj<length;jj++)
2115        {
2116                A1[jj] = (*intPointOfFacet)[jj];
[4199b3]2117                A2[jj] = -(*ivw)[jj];//TODO Or minus (*ivw)[ii] ??? NOTE minus
[be7ab3f]2118        }
2119        omFree(dstRing->wvhdl[0]);
2120        if(dstRing->wvhdl[1]!=NULL)
2121                omFree(dstRing->wvhdl[1]);
[4199b3]2122        dstRing->wvhdl[0]=(int*)A1;             
[be7ab3f]2123        dstRing->block1[0]=length;
[4199b3]2124        dstRing->wvhdl[1]=(int*)A2;             
[be7ab3f]2125        dstRing->block1[1]=length;
2126        rComplete(dstRing);
2127        rChangeCurrRing(dstRing);
[4199b3]2128        ideal dstRing_I;                       
[be7ab3f]2129        dstRing_I=idrCopyR(srcRing_HH,srcRing);
[4199b3]2130        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more       
[be7ab3f]2131        save=test;
2132        test|=Sy_bit(OPT_REDSB);
2133        test|=Sy_bit(OPT_REDTAIL);
2134        ideal tmpI;
2135        tmpI = dstRing_I;
2136#ifdef gfanp
[4199b3]2137//      timeval t_kStd_start, t_kStd_end;
[be7ab3f]2138        gettimeofday(&t_kStd_start,0);
2139#endif
[4199b3]2140//      if(gcone::hasHomInput==TRUE)
2141//              dstRing_I=kStd(tmpI,NULL,isHomog,NULL/*,gcone::hilbertFunction*/);
2142//      else
[be7ab3f]2143                dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
2144#ifdef gfanp
2145        gettimeofday(&t_kStd_end, 0);
2146        t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
2147#endif
2148        idDelete(&tmpI);
[4199b3]2149        idNorm(dstRing_I);                     
[be7ab3f]2150        idSkipZeroes(dstRing_I);
2151        test=save;
2152        /*End of step 3 - reduction*/
[4199b3]2153       
[be7ab3f]2154        f->setFlipGB(dstRing_I);
2155        f->flipRing=rCopy(dstRing);
2156        rDelete(tmpRing);
2157        rDelete(dstRing);
[4199b3]2158        //Now we should have dstRing with (a(),a(),dp,C)
2159        //This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
2160        //of cones in noRevS
[be7ab3f]2161        rChangeCurrRing(srcRing);
2162#ifdef gfanp
2163        gettimeofday(&end, 0);
2164        time_flip2 += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2165#endif
[4199b3]2166}//flip2
[be7ab3f]2167
[7c0ec6]2168/** \brief Compute initial ideal
2169 * Compute the initial ideal in_v(G) wrt a (possible) facet normal
2170 * used in gcone::getFacetNormal in order to preprocess possible facet normals
2171 * and in gcone::flip for obvious reasons.
2172*/
[f91fddc]2173/*inline*/ void gcone::computeInv(const ideal &gb, ideal &initialForm, const int64vec &fNormal)
[7c0ec6]2174{
[42660f]2175#ifdef gfanp
2176        timeval start, end;
2177        gettimeofday(&start, 0);
2178#endif
[7c0ec6]2179        for (int ii=0;ii<IDELEMS(gb);ii++)
2180        {
[42660f]2181                poly initialFormElement;
[871f01]2182                poly aktpoly = (poly)gb->m[ii];//Ptr, so don't pDelete(aktpoly)
[7c0ec6]2183                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
[4199b3]2184                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
[dcf8b4b]2185                initialFormElement=pHead(aktpoly);
[4199b3]2186//              int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
[7c0ec6]2187                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
2188                {
[590f813]2189                        int64vec *check = new int64vec(this->numVars);
[4199b3]2190                        aktpoly=pNext(aktpoly); //next term
[42660f]2191                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
[4199b3]2192                        pGetExpV(aktpoly,v);           
[f91fddc]2193                        /* Convert (int)v into (int64vec)check */
[4199b3]2194//                      bool notPar=FALSE;
[7c0ec6]2195                        for (int jj=0;jj<this->numVars;jj++)
2196                        {
2197                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
[4199b3]2198//                              register int64 foo=(fNormal)[jj];
2199//                              if( ( (*check)[jj]  == /*fNormal[jj]*/foo )
2200//                               || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) )
2201//                              {
2202//                                      notPar=TRUE;
2203//                                      break;
2204//                              }
[7c0ec6]2205                        }
[f91fddc]2206                        omFree(v);
[4199b3]2207                        if (isParallel(*check,fNormal))//Found a parallel vector. Add it
2208//                      if(notPar==FALSE)
[7c0ec6]2209                        {
[4199b3]2210                                initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));//pAdd = p_Add_q destroys args
[42660f]2211                        }
2212                        delete check;
[4199b3]2213                }//while
2214//              omFree(v);
[96e408]2215#ifndef NDEBUG
[4199b3]2216//              cout << "Initial Form=";                               
2217//              pWrite(initialFormElement[ii]);
2218//              cout << "---" << endl;
[7c0ec6]2219#endif
2220                /*Now initialFormElement must be added to (ideal)initialForm */
[dcf8b4b]2221                initialForm->m[ii]=pCopy(initialFormElement);
2222                pDelete(&initialFormElement);
[7c0ec6]2223                omFree(leadExpV);
[4199b3]2224        }//for
[42660f]2225#ifdef gfanp
2226        gettimeofday(&end, 0);
2227        time_computeInv += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2228#endif
[7c0ec6]2229}
2230
[1d9469]2231/** \brief Compute the remainder of a polynomial by a given ideal
2232 *
2233 * Compute \f$ f^{\mathcal{G}} \f$
2234 * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
2235 * However, since we are only interested in the remainder, there is no need to
2236 * compute the factors \f$ a_i \f$
2237 */
[4199b3]2238//NOTE: Should be replaced by kNF or kNF2
2239//NOTE: Done
2240//NOTE: removed with r12286
2241               
[1d9469]2242/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
[0d278cd]2243*/
[4199b3]2244//NOTE: use kNF or kNF2 instead of restOfDivision
[a23a297]2245inline ideal gcone::ffG(const ideal &H, const ideal &G)
[1d9469]2246{
2247        int size=IDELEMS(H);
2248        ideal res=idInit(size,1);
2249        for (int ii=0;ii<size;ii++)
2250        {
[4199b3]2251//              poly temp1;//=pInit();
2252//              poly temp2;//=pInit();
2253                poly temp3;//=pInit();//polys to temporarily store values for pSub
2254//              res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
2255//              temp1=pCopy(H->m[ii]);//TRY
2256//              temp2=pCopy(res->m[ii]);
2257                //NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
2258//              temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));//TRY
2259//              temp3=pSub(temp1, temp2);//TRY
2260                temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0)));//NOTRY
[6d0c34]2261                res->m[ii]=pCopy(temp3);
[4199b3]2262                //res->m[ii]=pSub(temp1,temp2); //buggy         
2263                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
2264//              pDelete(&temp1);//TRY
2265//              pDelete(&temp2);
[147167f]2266                pDelete(&temp3);
[42660f]2267        }
[1d9469]2268        return res;
2269}
[4199b3]2270       
[2e06300]2271/** \brief Preprocessing of inequlities
2272* Do some preprocessing on the matrix of inequalities
2273* 1) Replace several constraints on the pos. orthants by just one for each orthant
2274* 2) Remove duplicates of inequalities
2275* 3) Remove inequalities that arise as sums of other inequalities
2276*/
2277void gcone::preprocessInequalities(dd_MatrixPtr &ddineq)
2278{
[a7efc97]2279/*      int *posRowsArray=NULL;
[1d43d18]2280        int num_alloc=0;
2281        int num_elts=0;
[a7efc97]2282        int offset=0;*/
[4199b3]2283        //Remove zeroes (and strictly pos rows?)
[2e06300]2284        for(int ii=0;ii<ddineq->rowsize;ii++)
2285        {
[590f813]2286                int64vec *iv = new int64vec(this->numVars);
[4199b3]2287                int64vec *ivNull = new int64vec(this->numVars);//Needed for intvec64::compare(*int64vec)
[2e06300]2288                int posCtr=0;
2289                for(int jj=0;jj<this->numVars;jj++)
2290                {
2291                        (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
[4199b3]2292                        if((*iv)[jj]>0)//check for strictly pos rows
[a7efc97]2293                                posCtr++;
[4199b3]2294                        //Behold! This will delete the row for the standard simplex!
[2e06300]2295                }
[4199b3]2296//              if( (iv->compare(0)==0) || (posCtr==iv->length()) )
2297                if( (posCtr==iv->length()) || (iv->compare(ivNull)==0) )               
[2e06300]2298                {
2299                        dd_MatrixRowRemove(&ddineq,ii+1);
[4199b3]2300                        ii--;//Yes. This is on purpose
[2e06300]2301                }
2302                delete iv;
[1f3450]2303                delete ivNull;
[2e06300]2304        }
[4199b3]2305        //Remove duplicates of rows
2306//      posRowsArray=NULL;
2307//      num_alloc=0;
2308//      num_elts=0;
2309//      offset=0;
2310//      int num_newRows = ddineq->rowsize;
2311//      for(int ii=0;ii<ddineq->rowsize-1;ii++)
2312//      for(int ii=0;ii<num_newRows-1;ii++)
2313//      {
2314//              int64vec *iv = new int64vec(this->numVars);//1st vector to check against
2315//              for(int jj=0;jj<this->numVars;jj++)
2316//                      (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
2317//              for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++)
2318//              {
2319//                      int64vec *ivCheck = new int64vec(this->numVars);//Checked against iv
2320//                      for(int kk=0;kk<this->numVars;kk++)
2321//                              (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
2322//                      if (iv->compare(ivCheck)==0)
2323//                      {
2324// //                           cout << "=" << endl;
2325// //                           num_alloc++;
2326// //                           void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int)));
2327// //                           if(!tmp)
2328// //                           {
2329// //                                   WerrorS("Woah dude! Couldn't realloc memory\n");
2330// //                                   exit(-1);
2331// //                           }
2332// //                           posRowsArray = (int*)tmp;
2333// //                           posRowsArray[num_elts]=jj;
2334// //                           num_elts++;
2335//                              dd_MatrixRowRemove(&ddineq,jj+1);
2336//                              num_newRows = ddineq->rowsize;
2337//                      }
2338//                      delete ivCheck;
2339//              }
2340//              delete iv;
2341//      }
2342//      for(int ii=0;ii<num_elts;ii++)
2343//      {
2344//              dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
2345//              offset++;
2346//      }
2347//      free(posRowsArray);
2348        //Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/       
2349}//preprocessInequalities
2350       
[1d9469]2351/** \brief Compute a Groebner Basis
2352 *
2353 * Computes the Groebner basis and stores the result in gcone::gcBasis
2354 *\param ideal
2355 *\return void
2356 */
[4199b3]2357inline void gcone::getGB(const ideal &inputIdeal)               
[1d9469]2358{
2359        BITSET save=test;
2360        test|=Sy_bit(OPT_REDSB);
2361        test|=Sy_bit(OPT_REDTAIL);
2362        ideal gb;
2363        gb=kStd(inputIdeal,NULL,testHomog,NULL);
[bb503c7]2364        idNorm(gb);
[1d9469]2365        idSkipZeroes(gb);
[4199b3]2366        this->gcBasis=gb;       //write the GB into gcBasis
[1d9469]2367        test=save;
[4199b3]2368}//void getGB
2369               
[590f813]2370/** \brief Compute the negative of a given int64vec
[4199b3]2371        */             
[590f813]2372static int64vec* ivNeg(/*const*/int64vec *iv)
[4199b3]2373{       //Hm, switching to int64vec const int64vec does no longer work
2374        int64vec *res;// = new int64vec(iv->length());
[590f813]2375        res=iv64Copy(iv);
[4199b3]2376        *res *= (int)-1;                                               
[1d9469]2377        return res;
2378}
[e2fb7a]2379
[90adba8]2380
[1d9469]2381/** \brief Compute the dot product of two intvecs
2382*
2383*/
[4199b3]2384static int dotProduct(const int64vec &iva, const int64vec &ivb)                         
2385{                       
2386        int res=0;     
[1f637e]2387        for (int i=0;i<(currRing->N);i++)
[5e4025]2388        {
[4199b3]2389// #ifndef NDEBUG
2390//      (const_cast<int64vec*>(&iva))->show(1,0); (const_cast<int64vec*>(&ivb))->show(1,0);
2391// #endif
[5e4025]2392                res = res+(iva[i]*ivb[i]);
2393        }
2394        return res;
2395}
[1d9469]2396/** \brief Check whether two intvecs are parallel
2397 *
2398 * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
2399 */
[590f813]2400static bool isParallel(const int64vec &a,const int64vec &b)
[4199b3]2401{       
[1d9469]2402        bool res;
[f91fddc]2403        int lhs=dotProduct(a,b)*dotProduct(a,b);
[73f78d]2404        int rhs=dotProduct(a,a)*dotProduct(b,b);
2405        return res = (lhs==rhs)?TRUE:FALSE;
[e58dbb]2406}
[68eb0c]2407
[1d9469]2408/** \brief Compute an interior point of a given cone
[4199b3]2409 * Result will be written into int64vec iv.
[1d9469]2410 * Any rational point is automatically converted into an integer.
2411 */
[4199b3]2412void gcone::interiorPoint( dd_MatrixPtr &M, int64vec &iv) //no const &M here since we want to remove redundant rows
[1d9469]2413{
2414        dd_LPPtr lp,lpInt;
2415        dd_ErrorType err=dd_NoError;
2416        dd_LPSolverType solver=dd_DualSimplex;
2417        dd_LPSolutionPtr lpSol=NULL;
[4199b3]2418//      dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
2419//      dd_rowindex ddnewpos;
2420        dd_NumberType numb;     
2421        //M->representation=dd_Inequality;
2422                       
2423        //NOTE: Make this n-dimensional!
2424        //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
2425                                                                       
[2e06300]2426        /*NOTE: Leave the following line commented out!
2427        * Otherwise it will slow down computations a great deal
2428        * */
[4199b3]2429//      dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
2430        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
[2e06300]2431        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
2432        int jj=1;
2433        for (int ii=0;ii<this->numVars;ii++)
2434        {
2435                dd_set_si(posRestr->matrix[ii][jj],1);
[4199b3]2436                jj++;                                                   
[2e06300]2437        }
2438        dd_MatrixAppendTo(&M,posRestr);
2439        dd_FreeMatrix(posRestr);
[1d9469]2440        lp=dd_Matrix2LP(M, &err);
[a23a297]2441        if (err!=dd_NoError){WerrorS("Error during dd_Matrix2LP in gcone::interiorPoint");}
2442        if (lp==NULL){WerrorS("LP is NULL");}
[96e408]2443#ifndef NDEBUG
[4199b3]2444//                      dd_WriteLP(stdout,lp);
2445#endif 
2446                                       
[1d9469]2447        lpInt=dd_MakeLPforInteriorFinding(lp);
[a23a297]2448        if (err!=dd_NoError){WerrorS("Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint");}
[96e408]2449#ifndef NDEBUG
[4199b3]2450//                      dd_WriteLP(stdout,lpInt);
2451#endif                 
2452//      dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
[1d9469]2453        if (err!=dd_NoError)
2454        {
[a23a297]2455                WerrorS("Error during dd_FindRelativeInterior in gcone::interiorPoint");
[1d9469]2456                dd_WriteErrorMessages(stdout, err);
[4199b3]2457        }                       
2458        dd_LPSolve(lpInt,solver,&err);  //This will not result in a point from the relative interior
2459//      if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}                                                                       
[2e06300]2460        lpSol=dd_CopyLPSolution(lpInt);
[4199b3]2461//      if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
[96e408]2462#ifndef NDEBUG
[11a7dc]2463        printf("Interior point: ");
[1d9469]2464        for (int ii=1; ii<(lpSol->d)-1;ii++)
2465        {
2466                dd_WriteNumber(stdout,lpSol->sol[ii]);
2467        }
[11a7dc]2468        printf("\n");
[4199b3]2469#endif                 
2470        //NOTE The following strongly resembles parts of makeInt.
2471        //Maybe merge sometimes
[1d9469]2472        mpz_t kgV; mpz_init(kgV);
2473        mpz_set_str(kgV,"1",10);
2474        mpz_t den; mpz_init(den);
2475        mpz_t tmp; mpz_init(tmp);
2476        mpq_get_den(tmp,lpSol->sol[1]);
2477        for (int ii=1;ii<(lpSol->d)-1;ii++)
2478        {
2479                mpq_get_den(den,lpSol->sol[ii+1]);
2480                mpz_lcm(kgV,tmp,den);
2481                mpz_set(tmp, kgV);
2482        }
2483        mpq_t qkgV;
2484        mpq_init(qkgV);
2485        mpq_set_z(qkgV,kgV);
2486        for (int ii=1;ii<(lpSol->d)-1;ii++)
2487        {
2488                mpq_t product;
2489                mpq_init(product);
2490                mpq_mul(product,qkgV,lpSol->sol[ii]);
2491                iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
2492                mpq_clear(product);
2493        }
[96e408]2494#ifndef NDEBUG
[4199b3]2495//                      iv.show();
2496//                      cout << endl;
[ae5056]2497#endif
[1d9469]2498        mpq_clear(qkgV);
2499        mpz_clear(tmp);
2500        mpz_clear(den);
[4199b3]2501        mpz_clear(kgV);                 
2502                       
[1d9469]2503        dd_FreeLPSolution(lpSol);
2504        dd_FreeLPData(lpInt);
2505        dd_FreeLPData(lp);
[4199b3]2506//      set_free(ddlinset);
2507//      set_free(ddredrows);   
2508                       
2509}//void interiorPoint(dd_MatrixPtr const &M)
[2e06300]2510
[be7ab3f]2511/** Computes an interior point of a cone by taking two interior points a,b from two different facets
2512* and then computing b+(a-b)/2
[4199b3]2513* Of course this only works for flippable facets
2514* Two cases may occur:
[be7ab3f]2515* 1st: There are only two facets who share the only strictly positive ray
2516* 2nd: There are at least two facets which have a distinct positive ray
2517* In the former case we use linear algebra to determine an interior point,
2518* in the latter case we simply add up the two rays
2519*
2520* Way too bad! The case may occur that the cone is spanned by three rays, of which only two are strictly
2521* positive => these lie in a plane and thus their sum is not from relative interior.
2522* So let's just sum up all rays, find one strictly positive and shift the point along that ray
[4199b3]2523*
[be7ab3f]2524* Used by noRevS
[147167f]2525*NOTE no longer used nor maintained. MM Mar 9, 2010
[be7ab3f]2526*/
[4199b3]2527// void gcone::interiorPoint2()
2528// {//idPrint(this->gcBasis);
[96e408]2529// #ifndef NDEBUG
[4199b3]2530//      if(this->ivIntPt!=NULL)
2531//              WarnS("Interior point already exists - ovrewriting!");
2532// #endif
2533//      facet *f1 = this->facetPtr;
2534//      facet *f2 = NULL;
2535//      int64vec *intF1=NULL;
2536//      while(f1!=NULL)
2537//      {
2538//              if(f1->isFlippable)
2539//              {
2540//                      facet *f1Ray = f1->codim2Ptr;
2541//                      while(f1Ray!=NULL)
2542//                      {
2543//                              const int64vec *check=f1Ray->getRef2FacetNormal();
2544//                              if(iv64isStrictlyPositive(check))
2545//                              {
2546//                                      intF1=iv64Copy(check);
2547//                                      break;
2548//                              }                               
2549//                              f1Ray=f1Ray->next;
2550//                      }
2551//              }
2552//              if(intF1!=NULL)
2553//                      break;
2554//              f1=f1->next;
2555//      }
2556//      if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
2557//              f2=f1->next;
2558//      else
2559//              f2=this->facetPtr;
2560//      if(intF1==NULL && hasHomInput==TRUE)
2561//      {
2562//              intF1 = new int64vec(this->numVars);
2563//              for(int ii=0;ii<this->numVars;ii++)
2564//                      (*intF1)[ii]=1;
2565//      }
2566//      assert(f1); assert(f2);
2567//      int64vec *intF2=f2->getInteriorPoint();
2568//      mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
2569//      mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
2570//      mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
2571//      for(int ii=0;ii<this->numVars;ii++)
2572//      {
2573//              mpq_init(qPosRay[ii]);
2574//              mpq_init(qIntPt[ii]);
2575//              mpq_init(qPosIntPt[ii]);
2576//      }       
2577//      //Compute a+((b-a)/2) && Convert intF1 to mpq
2578//      for(int ii=0;ii<this->numVars;ii++)
2579//      {
2580//              mpq_t a,b;
2581//              mpq_init(a); mpq_init(b);
2582//              mpq_set_si(a,(*intF1)[ii],1);
2583//              mpq_set_si(b,(*intF2)[ii],1);
2584//              mpq_t diff;
2585//              mpq_init(diff);
2586//              mpq_sub(diff,b,a);      //diff=b-a
2587//              mpq_t quot;
2588//              mpq_init(quot);
2589//              mpq_div_2exp(quot,diff,1);      //quot=diff/2=(b-a)/2
2590//              mpq_clear(diff);
2591//              //Don't be clever and reuse diff here
2592//              mpq_t sum; mpq_init(sum);
2593//              mpq_add(sum,b,quot);    //sum=b+quot=a+(b-a)/2
2594//              mpq_set(qIntPt[ii],sum);
2595//              mpq_clear(sum);
2596//              mpq_clear(quot);
2597//              mpq_clear(a); mpq_clear(b);
2598//              //Now for intF1
2599//              mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
2600//      }
2601//      //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
2602//      while(TRUE)
2603//      {       
2604//              bool success=FALSE;
2605//              int posCtr=0;   
2606//              for(int ii=0;ii<this->numVars;ii++)
2607//              {
2608//                      mpq_t sum; mpq_init(sum);
2609//                      mpq_add(sum,qPosRay[ii],qIntPt[ii]);
2610//                      mpq_set(qPosIntPt[ii],sum);
2611//                      mpq_clear(sum);
2612//                      if(mpq_sgn(qPosIntPt[ii])==1)
2613//                              posCtr++;
2614//              }
2615//              if(posCtr==this->numVars)//qPosIntPt > 0
2616//                      break;
2617//              else
2618//              {
2619//                      mpq_t qTwo; mpq_init(qTwo);
2620//                      mpq_set_ui(qTwo,2,1);
2621//                      for(int jj=0;jj<this->numVars;jj++)
2622//                      {
2623//                              mpq_t tmp; mpq_init(tmp);
2624//                              mpq_mul(tmp,qPosRay[jj],qTwo);
2625//                              mpq_set( qPosRay[jj], tmp);
2626//                              mpq_clear(tmp);
2627//                      }
2628//                      mpq_clear(qTwo);
2629//              }
2630//      }//while
2631//      //Now qPosIntPt ought to be >0, so convert back to int :D
2632//      /*Compute lcm of the denominators*/
2633//      mpz_t *denom = new mpz_t[this->numVars];
2634//      mpz_t tmp,kgV;
2635//      mpz_init(tmp); mpz_init(kgV);
2636//      for (int ii=0;ii<this->numVars;ii++)
2637//      {
2638//              mpz_t z;
2639//              mpz_init(z);
2640//              mpq_get_den(z,qPosIntPt[ii]);
2641//              mpz_init(denom[ii]);
2642//              mpz_set( denom[ii], z);
2643//              mpz_clear(z);                           
2644//      }
2645//             
2646//      mpz_set(tmp,denom[0]);
2647//      for (int ii=0;ii<this->numVars;ii++)
2648//      {
2649//              mpz_lcm(kgV,tmp,denom[ii]);
2650//              mpz_set(tmp,kgV);                               
2651//      }
2652//      mpz_clear(tmp);
2653//      /*Multiply the nominators by kgV*/
2654//      mpq_t qkgV,res;
2655//      mpq_init(qkgV);
2656//      mpq_canonicalize(qkgV);         
2657//      mpq_init(res);
2658//      mpq_canonicalize(res);
2659//                             
2660//      mpq_set_num(qkgV,kgV);
2661//      int64vec *n=new int64vec(this->numVars);
2662//      for (int ii=0;ii<this->numVars;ii++)
2663//      {
2664//              mpq_canonicalize(qPosIntPt[ii]);
2665//              mpq_mul(res,qkgV,qPosIntPt[ii]);
2666//              (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
2667//      }
2668//      this->setIntPoint(n);
2669//      delete n;
2670//      delete [] qPosIntPt;
2671//      delete [] denom;
2672//      delete [] qPosRay;
2673//      delete [] qIntPt;
2674//      mpz_clear(kgV);
2675//      mpq_clear(qkgV); mpq_clear(res);
2676// }
2677       
[1d9469]2678/** \brief Copy a ring and add a weighted ordering in first place
[4199b3]2679 *
[1d9469]2680 */
[4199b3]2681ring gcone::rCopyAndAddWeight(const ring &r, int64vec *ivw)                             
[1d9469]2682{
2683        ring res=rCopy0(r);
2684        int jj;
[4199b3]2685                       
[1d9469]2686        omFree(res->order);
[590f813]2687        res->order =(int *)omAlloc0(4*sizeof(int/*64*/));
[1d9469]2688        omFree(res->block0);
[590f813]2689        res->block0=(int *)omAlloc0(4*sizeof(int/*64*/));
[1d9469]2690        omFree(res->block1);
[590f813]2691        res->block1=(int *)omAlloc0(4*sizeof(int/*64*/));
[1d9469]2692        omfree(res->wvhdl);
[590f813]2693        res->wvhdl =(int **)omAlloc0(4*sizeof(int/*64*/**));
[4199b3]2694                       
[590f813]2695        res->order[0]=ringorder_a/*64*/;
[1d9469]2696        res->block0[0]=1;
[590f813]2697        res->block1[0]=res->N;
[4199b3]2698        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                   
[1d9469]2699        res->block0[1]=1;
[590f813]2700        res->block1[1]=res->N;
[1d9469]2701        res->order[2]=ringorder_C;
[3449c9b]2702
[1d9469]2703        int length=ivw->length();
[590f813]2704        int/*64*/ *A=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
[1d9469]2705        for (jj=0;jj<length;jj++)
[4199b3]2706        {                               
[590f813]2707                A[jj]=(*ivw)[jj];
[dd6b1d]2708                if((*ivw)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::rCopyAndAddWeight!\n");
[4199b3]2709        }                       
[1d9469]2710        res->wvhdl[0]=(int *)A;
2711        res->block1[0]=length;
[4199b3]2712                       
[1d9469]2713        rComplete(res);
2714        return res;
[4199b3]2715}//rCopyAndAdd
2716               
2717ring gcone::rCopyAndAddWeight2(const ring &r,const int64vec *ivw, const int64vec *fNormal)                             
[be7ab3f]2718{
2719        ring res=rCopy0(r);
[4199b3]2720                       
[be7ab3f]2721        omFree(res->order);
[590f813]2722        res->order =(int *)omAlloc0(5*sizeof(int/*64*/));
[be7ab3f]2723        omFree(res->block0);
[590f813]2724        res->block0=(int *)omAlloc0(5*sizeof(int/*64*/));
[be7ab3f]2725        omFree(res->block1);
[590f813]2726        res->block1=(int *)omAlloc0(5*sizeof(int/*64*/));
[be7ab3f]2727        omfree(res->wvhdl);
[590f813]2728        res->wvhdl =(int **)omAlloc0(5*sizeof(int/*64*/**));
[4199b3]2729                       
[590f813]2730        res->order[0]=ringorder_a/*64*/;
[be7ab3f]2731        res->block0[0]=1;
2732        res->block1[0]=res->N;
[590f813]2733        res->order[1]=ringorder_a/*64*/;
[be7ab3f]2734        res->block0[1]=1;
2735        res->block1[1]=res->N;
[4199b3]2736       
[be7ab3f]2737        res->order[2]=ringorder_dp;
2738        res->block0[2]=1;
2739        res->block1[2]=res->N;
[4199b3]2740       
[be7ab3f]2741        res->order[3]=ringorder_C;
2742
2743        int length=ivw->length();
[590f813]2744        int/*64*/ *A1=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
2745        int/*64*/ *A2=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
[be7ab3f]2746        for (int jj=0;jj<length;jj++)
[4199b3]2747        {                               
[be7ab3f]2748                A1[jj]=(*ivw)[jj];
2749                A2[jj]=-(*fNormal)[jj];
[dd6b1d]2750                if((*ivw)[jj]>=INT_MAX || (*fNormal)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::rCopyAndAddWeight2!\n");
[4199b3]2751        }                       
[be7ab3f]2752        res->wvhdl[0]=(int *)A1;
2753        res->block1[0]=length;
2754        res->wvhdl[1]=(int *)A2;
2755        res->block1[1]=length;
2756        rComplete(res);
2757        return res;
2758}
[4199b3]2759               
2760//NOTE not needed anywhere
2761// ring rCopyAndChangeWeight(ring const &r, int64vec *ivw)
2762// {
2763//      ring res=rCopy0(currRing);
2764//      rComplete(res);
2765//      rSetWeightVec(res,(int64*)ivw);
2766//      //rChangeCurrRing(rnew);
2767//      return res;
2768// }
2769               
[1d9469]2770/** \brief Checks whether a given facet is a search facet
2771 * Determines whether a given facet of a cone is the search facet of a neighbouring cone
2772 * This is done in the following way:
2773 * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
2774 * that is first crossed during the generic walk.
2775 * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
[4199b3]2776 * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
[1d9469]2777*/
[4199b3]2778//removed with r12286
2779               
[1d9469]2780/** \brief Check for equality of two intvecs
2781 */
[590f813]2782static bool ivAreEqual(const int64vec &a, const int64vec &b)
[1d9469]2783{
2784        bool res=TRUE;
[1f637e]2785        for(int ii=0;ii<(currRing->N);ii++)
[1d9469]2786        {
2787                if(a[ii]!=b[ii])
[eaa679]2788                {
[1d9469]2789                        res=FALSE;
2790                        break;
[eaa679]2791                }
[1d9469]2792        }
2793        return res;
2794}
[4199b3]2795               
[1d9469]2796/** \brief The reverse search algorithm
2797 */
[4199b3]2798//removed with r12286
[72e467]2799/** \brief Compute the lineality/homogeneity space
2800* It is the kernel of the inequality matrix Ax=0
[276932]2801* As a result gcone::dd_LinealitySpace is set
[72e467]2802*/
2803dd_MatrixPtr gcone::computeLinealitySpace()
2804{
2805        dd_MatrixPtr res;
2806        dd_MatrixPtr ddineq;
[4199b3]2807        ddineq=dd_CopyMatrix(this->ddFacets);   
2808        //Add a row of 0s in 0th place
[72e467]2809        dd_MatrixPtr ddAppendRowOfZeroes=dd_CreateMatrix(1,this->numVars+1);
2810        dd_MatrixPtr ddFoo=dd_AppendMatrix(ddAppendRowOfZeroes,ddineq);
2811        dd_FreeMatrix(ddAppendRowOfZeroes);
2812        dd_FreeMatrix(ddineq);
2813        ddineq=dd_CopyMatrix(ddFoo);
2814        dd_FreeMatrix(ddFoo);
[4199b3]2815        //Cohen starts here
2816        int dimKer=0;//Cohen calls this r
2817        int m=ddineq->rowsize;//Rows
2818        int n=ddineq->colsize;//Cols
[72e467]2819        int c[m+1];
2820        int d[n+1];
2821        for(int ii=0;ii<m;ii++)
2822                c[ii]=0;
2823        for(int ii=0;ii<n;ii++)
[4199b3]2824                d[ii]=0;       
2825       
[72e467]2826        for(int k=1;k<n;k++)
2827        {
[4199b3]2828                //Let's find a j s.t. m[j][k]!=0 && c[j]=0             
2829                int condCtr=0;//Check each row for zeroness
[72e467]2830                for(int j=1;j<m;j++)
2831                {
2832                        if(mpq_sgn(ddineq->matrix[j][k])!=0 && c[j]==0)
2833                        {
2834                                mpq_t quot; mpq_init(quot);
2835                                mpq_t one; mpq_init(one); mpq_set_str(one,"-1",10);
2836                                mpq_t ratd; mpq_init(ratd);
2837                                if((int)mpq_get_d(ddineq->matrix[j][k])!=0)
2838                                        mpq_div(quot,one,ddineq->matrix[j][k]);
2839                                mpq_set(ratd,quot);
2840                                mpq_canonicalize(ratd);
[4199b3]2841               
[72e467]2842                                mpq_set_str(ddineq->matrix[j][k],"-1",10);
2843                                for(int ss=k+1;ss<n;ss++)
2844                                {
2845                                        mpq_t prod; mpq_init(prod);
[4199b3]2846                                        mpq_mul(prod, ratd, ddineq->matrix[j][ss]);                             
[72e467]2847                                        mpq_set(ddineq->matrix[j][ss],prod);
2848                                        mpq_canonicalize(ddineq->matrix[j][ss]);
2849                                        mpq_clear(prod);
[4199b3]2850                                }               
[72e467]2851                                for(int ii=1;ii<m;ii++)
2852                                {
2853                                        if(ii!=j)
2854                                        {
2855                                                mpq_set(ratd,ddineq->matrix[ii][k]);
[4199b3]2856                                                mpq_set_str(ddineq->matrix[ii][k],"0",10);                     
[1f3450]2857                                                for(int ss=k+1;ss<n;ss++)
[72e467]2858                                                {
2859                                                        mpq_t prod; mpq_init(prod);
2860                                                        mpq_mul(prod, ratd, ddineq->matrix[j][ss]);
2861                                                        mpq_t sum; mpq_init(sum);
2862                                                        mpq_add(sum, ddineq->matrix[ii][ss], prod);
2863                                                        mpq_set(ddineq->matrix[ii][ss], sum);
2864                                                        mpq_canonicalize(ddineq->matrix[ii][ss]);
2865                                                        mpq_clear(prod);
2866                                                        mpq_clear(sum);
2867                                                }
2868                                        }
2869                                }
[4199b3]2870                                c[j]=k;         
[72e467]2871                                d[k]=j;
[4199b3]2872                                mpq_clear(quot); mpq_clear(ratd); mpq_clear(one);       
[72e467]2873                        }
2874                        else
2875                                condCtr++;
[4199b3]2876                }                       
2877                if(condCtr==m-1)        //Why -1 ???
[72e467]2878                {
2879                        dimKer++;
2880                        d[k]=0;
[4199b3]2881//                      break;//goto _4;
2882                }//else{
[72e467]2883                /*Eliminate*/
[4199b3]2884        }//for(k
2885        /*Output kernel, i.e. set gcone::dd_LinealitySpace here*/       
2886//      gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
[4c4a80]2887        res = dd_CreateMatrix(dimKer,this->numVars+1);
[72e467]2888        int row=-1;
2889        for(int kk=1;kk<n;kk++)
2890        {
2891                if(d[kk]==0)
2892                {
2893                        row++;
2894                        for(int ii=1;ii<n;ii++)
2895                        {
2896                                if(d[ii]>0)
[4c4a80]2897                                        mpq_set(res->matrix[row][ii],ddineq->matrix[d[ii]][kk]);
[4199b3]2898                                else if(ii==kk)                         
[4c4a80]2899                                        mpq_set_str(res->matrix[row][ii],"1",10);
2900                                mpq_canonicalize(res->matrix[row][ii]);
[72e467]2901                        }
2902                }
2903        }
2904        dd_FreeMatrix(ddineq);
[4c4a80]2905        return(res);
[4199b3]2906        //Better safe than sorry:
2907        //NOTE dd_LinealitySpace contains RATIONAL ENTRIES
2908        //Therefore if you need integer ones: CALL gcone::makeInt(...) method
2909}       
[4c4a80]2910
2911
[1d9469]2912/** \brief The new method of Markwig and Jensen
2913 * Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
2914 * In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
[590f813]2915 * Keep a list of facets as a linked list containing an int64vec and an integer matrix.
[1d9469]2916 * Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
2917 * soon as we compute this facet again. Comparison of facets is done by...
2918 */
2919void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
[4199b3]2920{       
2921        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
[1d9469]2922        facet *SearchListAct;
2923        SearchListAct = NULL;
[4199b3]2924        //SearchListAct = SearchListRoot;                       
[1d9469]2925        gcone *gcAct;
2926        gcAct = &gcRoot;
[4199b3]2927        gcone *gcPtr;   //Pointer to end of linked list of cones
[1d9469]2928        gcPtr = &gcRoot;
[4199b3]2929        gcone *gcNext;  //Pointer to next cone to be visited
[1d9469]2930        gcNext = &gcRoot;
2931        gcone *gcHead;
[4199b3]2932        gcHead = &gcRoot;                       
[1d9469]2933        facet *fAct;
[4199b3]2934        fAct = gcAct->facetPtr;                         
[1d9469]2935        ring rAct;
2936        rAct = currRing;
[4199b3]2937        int UCNcounter=gcAct->getUCN(); 
[96e408]2938#ifndef NDEBUG
[11a7dc]2939        printf("NoRevs\n");
2940        printf("Facets are:\n");
[1d9469]2941        gcAct->showFacets();
[4199b3]2942#endif                 
[72e467]2943        /*Compute lineality space here
[276932]2944        * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/
[be7ab3f]2945        if(dd_LinealitySpace==NULL)
2946                dd_LinealitySpace = gcAct->computeLinealitySpace();
[4199b3]2947        /*End of lineality space computation*/         
2948        //gcAct->getCodim2Normals(*gcAct);
[be7ab3f]2949        if(fAct->codim2Ptr==NULL)
[4199b3]2950                gcAct->getExtremalRays(*gcAct);         
[963eeb]2951        /* Make a copy of the facet list of first cone
2952           Since the operations getCodim2Normals and normalize affect the facets
[4199b3]2953           we must not memcpy them before these ops! */ 
2954        /*facet *codim2Act; codim2Act = NULL;                   
[276932]2955        facet *sl2Root;
[4199b3]2956        facet *sl2Act;*/                       
[1d9469]2957        for(int ii=0;ii<this->numFacets;ii++)
2958        {
[4199b3]2959                //only copy flippable facets! NOTE: We do not need flipRing or any such information.
[1d9469]2960                if(fAct->isFlippable==TRUE)
2961                {
[73809b]2962                        /*Using shallow copy here*/
2963#ifdef SHALLOW
[4199b3]2964                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
[73809b]2965                        {
[64b6c6]2966                                if(SearchListRoot!=NULL) delete(SearchListRoot);
[73809b]2967                                SearchListRoot = fAct->shallowCopy(*fAct);
[4199b3]2968                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
[73809b]2969                        }
2970                        else
[4199b3]2971                        {                       
[73809b]2972                                SearchListAct->next = fAct->shallowCopy(*fAct);
[4199b3]2973                                SearchListAct = SearchListAct->next;                                           
[73809b]2974                        }
2975                        fAct=fAct->next;
2976#else
2977                        /*End of shallow copy*/
[590f813]2978                        int64vec *fNormal;
[1d9469]2979                        fNormal = fAct->getFacetNormal();
[4199b3]2980                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
2981                        {                                               
2982                                SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
[1d9469]2983                        }
2984                        else
[4199b3]2985                        {                       
[1d9469]2986                                SearchListAct->next = new facet();
[4199b3]2987                                SearchListAct = SearchListAct->next;                                           
[1d9469]2988                        }
2989                        SearchListAct->setFacetNormal(fNormal);
2990                        SearchListAct->setUCN(this->getUCN());
2991                        SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
2992                        SearchListAct->isFlippable=TRUE;
[4199b3]2993                        //Copy int point as well
[64b6c6]2994                        int64vec *ivIntPt;
[be7ab3f]2995                        ivIntPt = fAct->getInteriorPoint();
2996                        SearchListAct->setInteriorPoint(ivIntPt);
2997                        delete(ivIntPt);
[4199b3]2998                        //Copy codim2-facets
[276932]2999                        facet *codim2Act;
[4199b3]3000                        codim2Act = NULL;                       
[276932]3001                        facet *sl2Root;
[4199b3]3002                        facet *sl2Act;                 
[1d9469]3003                        codim2Act=fAct->codim2Ptr;
3004                        SearchListAct->codim2Ptr = new facet(2);
3005                        sl2Root = SearchListAct->codim2Ptr;
[4199b3]3006                        sl2Act = sl2Root;                       
[1d9469]3007                        for(int jj=0;jj<fAct->numCodim2Facets;jj++)
[4199b3]3008//                      for(int jj=0;jj<fAct->numRays-1;jj++)
[2b3c62f]3009                        {
[590f813]3010                                int64vec *f2Normal;
[1d9469]3011                                f2Normal = codim2Act->getFacetNormal();
3012                                if(jj==0)
[4199b3]3013                                {                                               
[1d9469]3014                                        sl2Act = sl2Root;
3015                                        sl2Act->setFacetNormal(f2Normal);
[2b3c62f]3016                                }
[1d9469]3017                                else
3018                                {
3019                                        sl2Act->next = new facet(2);
3020                                        sl2Act = sl2Act->next;
3021                                        sl2Act->setFacetNormal(f2Normal);
[dcf8b4b]3022                                }
[4199b3]3023                                delete f2Normal;                               
[1d9469]3024                                codim2Act = codim2Act->next;
[2b3c62f]3025                        }
[1d9469]3026                        fAct = fAct->next;
[73809b]3027                        delete fNormal;
3028#endif
[4199b3]3029                }//if(fAct->isFlippable==TRUE)                 
[1d9469]3030                else {fAct = fAct->next;}
[4199b3]3031        }//End of copying facets into SLA
3032       
3033        SearchListAct = SearchListRoot; //Set to beginning of list
[1d9469]3034        /*Make SearchList doubly linked*/
[73809b]3035#ifndef NDEBUG
3036  #if SIZEOF_LONG==8
3037        while(SearchListAct!=(facet*)0xfefefefefefefefe && SearchListAct!=NULL)
3038        {
3039                if(SearchListAct->next!=(facet*)0xfefefefefefefefe && SearchListAct->next!=NULL){
3040  #elif SIZEOF_LONG!=8
3041        while(SearchListAct!=(facet*)0xfefefefe)
3042        {
3043                if(SearchListAct->next!=(facet*)0xfefefefe){
3044  #endif
3045#else
[1d9469]3046        while(SearchListAct!=NULL)
3047        {
[73809b]3048                if(SearchListAct->next!=NULL){
3049#endif
[4199b3]3050                        SearchListAct->next->prev = SearchListAct;                                     
[1d9469]3051                }
3052                SearchListAct = SearchListAct->next;
3053        }
[4199b3]3054        SearchListAct = SearchListRoot; //Set to beginning of List
3055       
3056//      fAct = gcAct->facetPtr;//???
3057        gcAct->writeConeToFile(*gcAct);                 
[1d9469]3058        /*End of initialisation*/
[4199b3]3059       
[1d9469]3060        fAct = SearchListAct;
[963eeb]3061        /*2nd step
[276932]3062          Choose a facet from SearchList, flip it and forget the previous cone
3063          We always choose the first facet from SearchList as facet to be flipped
[4199b3]3064        */     
3065        while( (SearchListAct!=NULL))//&& counter<490)
3066        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
3067                fAct = SearchListAct;           
[1d37196]3068                while(fAct!=NULL)
[4199b3]3069//              while( (fAct->getUCN() == fAct->next->getUCN()) )               
3070                {       //Since SLA should only contain flippables there should be no need to check for that                   
[be7ab3f]3071                        gcAct->flip2(gcAct->gcBasis,fAct);
[4199b3]3072                        //NOTE rCopy needed?
[1d9469]3073                        ring rTmp=rCopy(fAct->flipRing);
[63a2f8]3074                        rComplete(rTmp);
[1d9469]3075                        rChangeCurrRing(rTmp);
[b9a447]3076                        gcone *gcTmp = new gcone(*gcAct,*fAct);//copy constructor!
[bec3b0]3077                        /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
3078                         * Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
3079                         * NOTE: Can this cause trouble with the destructor?
3080                         * Answer: Yes it bloody well does!
3081                         * However I'd like to delete this data here, since if we have a cone with many many facets it
3082                         * should pay to get rid of the flipGB as soon as possible.
3083                         * Destructor must be equipped with necessary checks.
[963eeb]3084                        */
[6d0c34]3085                        idDelete((ideal *)&fAct->flipGB);
[4199b3]3086                        rDelete(fAct->flipRing);                       
[64b6c6]3087                        gcTmp->getConeNormals(gcTmp->gcBasis/*, FALSE*/);
[4199b3]3088//                      gcTmp->getCodim2Normals(*gcTmp);
[72e467]3089                        gcTmp->getExtremalRays(*gcTmp);
[4199b3]3090                        //NOTE Order rays lex here
3091                        gcTmp->orderRays();                     
3092//                      //NOTE If flip2 is used we need to get an interior point of gcTmp
3093//                      // and replace gcTmp->baseRing with an appropriate ring with only
3094//                      // one weight
3095//                      gcTmp->interiorPoint2();
[be7ab3f]3096                        gcTmp->replaceDouble_ringorder_a_ByASingleOne();
[4199b3]3097                        //gcTmp->ddFacets should not be needed anymore, so
[6d0c34]3098                        dd_FreeMatrix(gcTmp->ddFacets);
[96e408]3099#ifndef NDEBUG
[4199b3]3100//                      gcTmp->showFacets(1);
[195b03]3101#endif
[f07b38c]3102                        /*add facets to SLA here*/
[73809b]3103#ifdef SHALLOW
[96e408]3104  #ifndef NDEBUG
[64b6c6]3105                        printf("fActUCN before enq2: %i\n",fAct->getUCN());
3106  #endif
3107                        facet *tmp;
3108                        tmp=gcTmp->enqueue2(SearchListRoot);
[96e408]3109  #ifndef NDEBUG
[64b6c6]3110                        printf("\nheadUCN=%i\n",tmp->getUCN());
3111                        printf("fActUCN after enq2: %i\n",fAct->getUCN());
3112  #endif
[63a2f8]3113                        SearchListRoot=tmp;
[4199b3]3114                        //SearchListRoot=gcTmp->enqueue2/*NewFacets*/(SearchListRoot);
3115#else
[dcf8b4b]3116                        SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
[4199b3]3117#endif //ifdef SHALLOW
3118//                      gcTmp->writeConeToFile(*gcTmp);
[d2cd515]3119                        if(gfanHeuristic==1)
[bec3b0]3120                        {
[d19585]3121                                gcTmp->writeConeToFile(*gcTmp);
[4199b3]3122                                idDelete((ideal*)&gcTmp->gcBasis);//Whonder why?
[1f3450]3123                                rDelete(gcTmp->baseRing);
[4199b3]3124                        }                       
[96e408]3125#ifndef NDEBUG
[1f3450]3126                        if(SearchListRoot!=NULL)
[11a7dc]3127                                showSLA(*SearchListRoot);
[4199b3]3128#endif                 
[1d9469]3129                        rChangeCurrRing(gcAct->baseRing);
[42660f]3130                        rDelete(rTmp);
[4199b3]3131                        //doubly linked for easier removal
[c5940e]3132                        gcTmp->prev = gcPtr;
[1d9469]3133                        gcPtr->next=gcTmp;
3134                        gcPtr=gcPtr->next;
[4199b3]3135                        //Cleverly disguised exit condition follows
[1d9469]3136                        if(fAct->getUCN() == fAct->next->getUCN())
3137                        {
[63a2f8]3138                                printf("Switching UCN from %i to %i\n",fAct->getUCN(),fAct->next->getUCN());
[4199b3]3139                                fAct=fAct->next;                               
[1d9469]3140                        }
3141                        else
[63a2f8]3142                        {
[4199b3]3143                                //rDelete(gcAct->baseRing);
3144//                              printf("break\n");
[1d9469]3145                                break;
[63a2f8]3146                        }
[4199b3]3147//                      fAct=fAct->next;
3148                }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );         
3149                //Search for cone with smallest UCN
[1f3450]3150#ifndef NDEBUG
[4199b3]3151  #if SIZEOF_LONG==8    //64 bit
[1f3450]3152                while(gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && SearchListRoot!=NULL)
[276932]3153  #elif SIZEOF_LONG == 4
[1f3450]3154                while(gcNext!=(gcone * const)0xfbfbfbfb && SearchListRoot!=NULL)
3155  #endif
3156#endif
3157#ifdef NDEBUG
[4199b3]3158                while(gcNext!=NULL && SearchListRoot!=NULL)     
[1f3450]3159#endif
[4199b3]3160                {                       
[1d37196]3161                        if( gcNext->getUCN() == SearchListRoot->getUCN() )
[1f3450]3162                        {
[d2cd515]3163                                if(gfanHeuristic==0)
[bec3b0]3164                                {
3165                                        gcAct = gcNext;
[4199b3]3166                                        //Seems better to not to use rCopy here
3167//                                      rAct=rCopy(gcAct->baseRing);
[1d37196]3168                                        rAct=gcAct->baseRing;
[bec3b0]3169                                        rComplete(rAct);
[4199b3]3170                                        rChangeCurrRing(rAct);                                         
[bec3b0]3171                                        break;
3172                                }
[d2cd515]3173                                else if(gfanHeuristic==1)
[bec3b0]3174                                {
[2e06300]3175                                        gcone *gcDel;
[4199b3]3176                                        gcDel = gcAct;                                 
[bec3b0]3177                                        gcAct = gcNext;
[4199b3]3178                                        //Read st00f from file &
3179                                        //implant the GB into gcAct
[f07b38c]3180                                        readConeFromFile(gcAct->getUCN(), gcAct);
[4199b3]3181                                        //Kill the baseRing but ONLY if it is not the ring the computation started in!
3182//                                      if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
3183//                                              rDelete(gcDel->baseRing);
3184//                                      rAct=rCopy(gcAct->baseRing);
3185                                        /*The ring change occurs in readConeFromFile, so as to
[72e467]3186                                        assure that gcAct->gcBasis belongs to the right ring*/
[1d37196]3187                                        rAct=gcAct->baseRing;
[e98068]3188                                        rComplete(rAct);
3189                                        rChangeCurrRing(rAct);
[bec3b0]3190                                        break;
[4199b3]3191                                }                               
[276932]3192                        }
[63a2f8]3193                        else if(gcNext->getUCN() < SearchListRoot->getUCN() )
3194                        {
[4199b3]3195                                idDelete( (ideal*)&gcNext->gcBasis );                           
3196//                              rDelete(gcNext->baseRing);//TODO Why does this crash?
[63a2f8]3197                        }
[276932]3198                        /*else
3199                        {
3200                                if(gfanHeuristic==1)
3201                                {
3202                                        gcone *gcDel;
3203                                        gcDel = gcNext;
3204                                        if(gcDel->getUCN()!=1)
3205                                                rDelete(gcDel->baseRing);
[4199b3]3206                                }                                       
3207                        }*/                     
[bec3b0]3208                        gcNext = gcNext->next;
[1d9469]3209                }
3210                UCNcounter++;
[4199b3]3211                SearchListAct = SearchListRoot;         
[1d9469]3212        }
[11a7dc]3213        printf("\nFound %i cones - terminating\n", counter);
[4199b3]3214}//void noRevS(gcone &gc)       
3215               
3216               
[1d9469]3217/** \brief Make a set of rational vectors into integers
3218 *
[73809b]3219 * We compute the lcm of the denominators and multiply with this to get integer values.
[147167f]3220 * If the gcd of the nominators > 1 we divide by the gcd => primitive vector.
[590f813]3221 * Expects a new int64vec as 3rd parameter
3222 * \param dd_MatrixPtr,int64vec
[1d9469]3223 */
[590f813]3224void gcone::makeInt(const dd_MatrixPtr &M, const int line, int64vec &n)
[4199b3]3225{                       
[dcf8b4b]3226        mpz_t *denom = new mpz_t[this->numVars];
[1d9469]3227        for(int ii=0;ii<this->numVars;ii++)
3228        {
3229                mpz_init_set_str(denom[ii],"0",10);
3230        }
[4199b3]3231               
[1d9469]3232        mpz_t kgV,tmp;
3233        mpz_init(kgV);
3234        mpz_init(tmp);
[4199b3]3235                       
[1d9469]3236        for (int ii=0;ii<(M->colsize)-1;ii++)
3237        {
3238                mpz_t z;
3239                mpz_init(z);
3240                mpq_get_den(z,M->matrix[line-1][ii+1]);
3241                mpz_set( denom[ii], z);
[4199b3]3242                mpz_clear(z);                           
[1d9469]3243        }
[4199b3]3244                                               
[1d9469]3245        /*Compute lcm of the denominators*/
3246        mpz_set(tmp,denom[0]);
3247        for (int ii=0;ii<(M->colsize)-1;ii++)
3248        {
3249                mpz_lcm(kgV,tmp,denom[ii]);
[4199b3]3250                mpz_set(tmp,kgV);                               
[1d9469]3251        }
[4199b3]3252        mpz_clear(tmp); 
[1d9469]3253        /*Multiply the nominators by kgV*/
3254        mpq_t qkgV,res;
3255        mpq_init(qkgV);
3256        mpq_set_str(qkgV,"1",10);
3257        mpq_canonicalize(qkgV);
[4199b3]3258                       
[1d9469]3259        mpq_init(res);
3260        mpq_set_str(res,"1",10);
3261        mpq_canonicalize(res);
[4199b3]3262                       
[1d9469]3263        mpq_set_num(qkgV,kgV);
[4199b3]3264                       
3265//      mpq_canonicalize(qkgV);
3266//      int ggT=1;
[1d9469]3267        for (int ii=0;ii<(M->colsize)-1;ii++)
3268        {
3269                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
[590f813]3270                n[ii]=(int64)mpz_get_d(mpq_numref(res));
[4199b3]3271//              ggT=int64gcd(ggT,n[ii]);
[73809b]3272        }
[590f813]3273        int64 ggT=n[0];
[e6ddbd]3274        for(int ii=0;ii<this->numVars;ii++)
[4199b3]3275                ggT=int64gcd(ggT,n[ii]);       
3276        //Normalisation
[73809b]3277        if(ggT>1)
3278        {
3279                for(int ii=0;ii<this->numVars;ii++)
3280                        n[ii] /= ggT;
[1d9469]3281        }
[dcf8b4b]3282        delete [] denom;
[1d9469]3283        mpz_clear(kgV);
[4199b3]3284        mpq_clear(qkgV); mpq_clear(res);                       
3285                       
[1d9469]3286}
3287/**
[4199b3]3288 * For all codim-2-facets we compute the gcd of the components of the facet normal and
[1d9469]3289 * divide it out. Thus we get a normalized representation of each
3290 * (codim-2)-facet normal, i.e. a primitive vector
[72e467]3291 * Actually we now also normalize the facet normals.
[1d9469]3292 */
[4199b3]3293// void gcone::normalize()
3294// {
3295//      int *ggT = new int;
3296//              *ggT=1;
3297//      facet *fAct;
3298//      facet *codim2Act;
3299//      fAct = this->facetPtr;
3300//      codim2Act = fAct->codim2Ptr;
3301//      while(fAct!=NULL)
3302//      {
3303//              int64vec *fNormal;
3304//              fNormal = fAct->getFacetNormal();
3305//              int *ggT = new int;
3306//              *ggT=1;
3307//              for(int ii=0;ii<this->numVars;ii++)
3308//              {
3309//                      *ggT=intgcd((*ggT),(*fNormal)[ii]);
3310//              }
3311//              if(*ggT>1)//We only need to do this if the ggT is non-trivial
3312//              {
3313// //                   int64vec *fCopy = fAct->getFacetNormal();
3314//                      for(int ii=0;ii<this->numVars;ii++)
3315//                              (*fNormal)[ii] = ((*fNormal)[ii])/(*ggT);
3316//                      fAct->setFacetNormal(fNormal);
3317//              }               
3318//              delete fNormal;
3319//              delete ggT;
3320//              /*And now for the codim2*/
3321//              while(codim2Act!=NULL)
3322//              {                               
3323//                      int64vec *n;
3324//                      n=codim2Act->getFacetNormal();
3325//                      int *ggT=new int;
3326//                      *ggT=1;
3327//                      for(int ii=0;ii<this->numVars;ii++)
3328//                      {
3329//                              *ggT = intgcd((*ggT),(*n)[ii]);
3330//                      }
3331//                      if(*ggT>1)
3332//                      {
3333//                              for(int ii=0;ii<this->numVars;ii++)
3334//                              {
3335//                                      (*n)[ii] = ((*n)[ii])/(*ggT);
3336//                              }
3337//                              codim2Act->setFacetNormal(n);
3338//                      }
3339//                      codim2Act = codim2Act->next;
3340//                      delete n;
3341//                      delete ggT;
3342//              }
3343//              fAct = fAct->next;
3344//      }
3345// }
3346
3347/** \brief Enqueue new facets into the searchlist
[1d9469]3348 * The searchlist (SLA for short) is implemented as a FIFO
3349 * Checks are done as follows:
3350 * 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
[4199b3]3351 * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
[1d9469]3352 *      facet from SLA and do not add fAct.
3353 * It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
3354 * can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA.
3355 * Takes ptr to search list root
3356 * Returns a pointer to new first element of Searchlist
3357 */
[1d37196]3358facet * gcone::enqueueNewFacets(facet *f)
[1d9469]3359{
[42660f]3360#ifdef gfanp
3361        timeval start, end;
3362        gettimeofday(&start, 0);
3363#endif
[1d9469]3364        facet *slHead;
[1d37196]3365        slHead = f;
[4199b3]3366        facet *slAct;   //called with f=SearchListRoot
[1d37196]3367        slAct = f;
[4199b3]3368        facet *slEnd;   //Pointer to end of SLA
[1d37196]3369        slEnd = f;
[4199b3]3370//      facet *slEndStatic;     //marks the end before a new facet is added             
[1d9469]3371        facet *fAct;
3372        fAct = this->facetPtr;
3373        facet *codim2Act;
3374        codim2Act = this->facetPtr->codim2Ptr;
3375        facet *sl2Act;
[1d37196]3376        sl2Act = f->codim2Ptr;
[1d9469]3377        /** Pointer to a facet that will be deleted */
[1d37196]3378        volatile facet *deleteMarker;
[1d9469]3379        deleteMarker = NULL;
[4199b3]3380                       
[1d9469]3381        /** \brief  Flag to mark a facet that might be added
3382         * The following case may occur:
3383         * A facet fAct is found to be parallel but not equal to the current facet slAct from SLA.
3384         * This does however not mean that there does not exist a facet behind the current slAct that is actually equal
3385         * to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
[4199b3]3386         * FALSE and the according slAct is deleted.
[1d9469]3387         * If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
3388         */
[64b6c6]3389
[4199b3]3390        /**A facet was removed, lengthOfSearchlist-- thus we must not rely on
[1d9469]3391         * if(notParallelCtr==lengthOfSearchList) but rather
3392         * if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
3393         */
3394        volatile bool removalOccured=FALSE;
3395        while(slEnd->next!=NULL)
3396        {
3397                slEnd=slEnd->next;
3398        }
[4199b3]3399        /*1st step: compare facetNormals*/                     
[1d9469]3400        while(fAct!=NULL)
[4199b3]3401        {                                               
[1d9469]3402                if(fAct->isFlippable==TRUE)
3403                {
[590f813]3404                        int64vec *fNormal=NULL;
[4199b3]3405                        fNormal=fAct->getFacetNormal();                 
[1d9469]3406                        slAct = slHead;
[4199b3]3407                        /*If slAct==NULL and fAct!=NULL
[1d9469]3408                        we just copy all remaining facets into SLA*/
3409                        if(slAct==NULL)
3410                        {
3411                                facet *fCopy;
3412                                fCopy = fAct;
3413                                while(fCopy!=NULL)
3414                                {
[4199b3]3415                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
[398e35]3416                                        {
[276932]3417                                                if(slAct==NULL)
3418                                                {
[4199b3]3419                                                        slAct = new facet(*fCopy/*,TRUE*/);//copy constructor
3420                                                        slHead = slAct;                                                         
[276932]3421                                                }
3422                                                else
3423                                                {
[73809b]3424                                                        slAct->next = new facet(*fCopy/*,TRUE*/);
[276932]3425                                                        slAct = slAct->next;
3426                                                }
[398e35]3427                                        }
[1d9469]3428                                        fCopy = fCopy->next;
3429                                }
[4199b3]3430                                break;//Where does this lead to?
[1d9469]3431                        }
[bec3b0]3432                        /*End of dumping into SLA*/
[1d9469]3433                        while(slAct!=NULL)
3434                        {
[590f813]3435                                int64vec *slNormal=NULL;
[1d9469]3436                                removalOccured=FALSE;
3437                                slNormal = slAct->getFacetNormal();
[96e408]3438#ifndef NDEBUG
[11a7dc]3439                                printf("Checking facet (");fNormal->show(1,1);printf(") against (");slNormal->show(1,1);printf(")\n");
[4199b3]3440#endif                         
3441//                              if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
3442//                                      exit(-1);
[73809b]3443                                if(areEqual2(fAct,slAct))
[4199b3]3444                                {                                       
[1d37196]3445                                        deleteMarker = slAct;
[f5a3167]3446                                        if(slAct==slHead)
[4199b3]3447                                        {                                               
3448                                                slHead = slAct->next;                                           
[f5a3167]3449                                                if(slHead!=NULL)
[1f3450]3450                                                        slHead->prev = NULL;
[f5a3167]3451                                        }
3452                                        else if (slAct==slEnd)
[1d9469]3453                                        {
[f5a3167]3454                                                slEnd=slEnd->prev;
[1f3450]3455                                                slEnd->next = NULL;
[4199b3]3456                                        }                                                               
[f5a3167]3457                                        else
[1d9469]3458                                        {
[f5a3167]3459                                                slAct->prev->next = slAct->next;
[1f3450]3460                                                slAct->next->prev = slAct->prev;
[f5a3167]3461                                        }
3462                                        removalOccured=TRUE;
[72e467]3463                                        gcone::lengthOfSearchList--;
[1d37196]3464                                        if(deleteMarker!=NULL)
3465                                        {
[4199b3]3466//                                              delete deleteMarker;
3467//                                              deleteMarker=NULL;
[1d37196]3468                                        }
[96e408]3469#ifndef NDEBUG
[11a7dc]3470                                        printf("Removing (");fNormal->show(1,1);printf(") from list\n");
[dcf8b4b]3471#endif
3472                                        delete slNormal;
[4199b3]3473                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
[1f3450]3474                                }
[1d9469]3475                                slAct = slAct->next;
[bec3b0]3476                                /* NOTE The following lines must not be here but rather called inside the if-block above.
[1d9469]3477                                Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now,
3478                                (not nice!) but since they are in seperate branches of the if-statement there should not
[bec3b0]3479                                be a way it gets called twice thus ommiting one facet:
[f5a3167]3480                                slAct = slAct->next;*/
3481                                if(deleteMarker!=NULL)
[4199b3]3482                                {                                               
3483//                                      delete deleteMarker;
3484//                                      deleteMarker=NULL;
[f5a3167]3485                                }
[dcf8b4b]3486                                delete slNormal;
[4199b3]3487                                                //if slAct was marked as to be deleted, delete it here!
3488                        }//while(slAct!=NULL)
[f5a3167]3489                        if(removalOccured==FALSE)
[1d9469]3490                        {
[96e408]3491#ifndef NDEBUG
[4199b3]3492//                              cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl;
[dcf8b4b]3493#endif
[4199b3]3494                                //Add fAct to SLA
[1d9469]3495                                facet *marker;
3496                                marker = slEnd;
3497                                facet *f2Act;
3498                                f2Act = fAct->codim2Ptr;
[73809b]3499
[1d9469]3500                                slEnd->next = new facet();
[4199b3]3501                                slEnd = slEnd->next;//Update slEnd
[1d9469]3502                                facet *slEndCodim2Root;
3503                                facet *slEndCodim2Act;
3504                                slEndCodim2Root = slEnd->codim2Ptr;
3505                                slEndCodim2Act = slEndCodim2Root;
[4199b3]3506                                                               
[1d9469]3507                                slEnd->setUCN(this->getUCN());
3508                                slEnd->isFlippable = TRUE;
3509                                slEnd->setFacetNormal(fNormal);
[4199b3]3510                                //NOTE Add interior point here
3511                                //This is ugly but needed for flip2
3512                                //Better: have slEnd->interiorPoint point to fAct->interiorPoint
3513                                //NOTE Only reference -> c.f. copy constructor
3514                                //slEnd->setInteriorPoint(fAct->getInteriorPoint());
[be7ab3f]3515                                slEnd->interiorPoint = fAct->interiorPoint;
[1d9469]3516                                slEnd->prev = marker;
[4199b3]3517                                //Copy codim2-facets
3518                                //int64vec *f2Normal=new int64vec(this->numVars);
[1d9469]3519                                while(f2Act!=NULL)
[dd2b206]3520                                {
[590f813]3521                                        int64vec *f2Normal;
[1d9469]3522                                        f2Normal=f2Act->getFacetNormal();
3523                                        if(slEndCodim2Root==NULL)
3524                                        {
3525                                                slEndCodim2Root = new facet(2);
[4199b3]3526                                                slEnd->codim2Ptr = slEndCodim2Root;                     
[1d9469]3527                                                slEndCodim2Root->setFacetNormal(f2Normal);
3528                                                slEndCodim2Act = slEndCodim2Root;
3529                                        }
[4199b3]3530                                        else                                   
[1d9469]3531                                        {
3532                                                slEndCodim2Act->next = new facet(2);
3533                                                slEndCodim2Act = slEndCodim2Act->next;
3534                                                slEndCodim2Act->setFacetNormal(f2Normal);
3535                                        }
3536                                        f2Act = f2Act->next;
[dcf8b4b]3537                                        delete f2Normal;
[dd2b206]3538                                }
[4199b3]3539                                gcone::lengthOfSearchList++;                                                   
3540                        }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
[1d9469]3541                        fAct = fAct->next;
[dcf8b4b]3542                        delete fNormal;
[4199b3]3543//                      delete slNormal;
3544                }//if(fAct->isFlippable==TRUE)
[1d9469]3545                else
3546                {
3547                        fAct = fAct->next;
3548                }
[72e467]3549                if(gcone::maxSize<gcone::lengthOfSearchList)
3550                        gcone::maxSize= gcone::lengthOfSearchList;
[4199b3]3551        }//while(fAct!=NULL)
[42660f]3552#ifdef gfanp
3553        gettimeofday(&end, 0);
3554        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
[4199b3]3555#endif                                         
[1d9469]3556        return slHead;
[4199b3]3557}//addC2N
[73809b]3558
[64b6c6]3559/** Enqueuing using shallow copies*/
[73809b]3560facet * gcone::enqueue2(facet *f)
3561{
3562        assert(f!=NULL);
3563#ifdef gfanp
3564        timeval start, end;
3565        gettimeofday(&start, 0);
3566#endif
3567        facet *slHead;
3568        slHead = f;
[4199b3]3569        facet *slAct;   //called with f=SearchListRoot
[73809b]3570        slAct = f;
[4199b3]3571        static facet *slEnd;    //Pointer to end of SLA
[147167f]3572        if(slEnd==NULL)
3573                slEnd = f;
[4199b3]3574       
[73809b]3575        facet *fAct;
[4199b3]3576        fAct = this->facetPtr;//New facets to compare
[73809b]3577        facet *codim2Act;
3578        codim2Act = this->facetPtr->codim2Ptr;
3579        volatile bool removalOccured=FALSE;
3580        while(slEnd->next!=NULL)
3581        {
3582                slEnd=slEnd->next;
[4199b3]3583        }       
[73809b]3584        while(fAct!=NULL)
3585        {
3586                if(fAct->isFlippable)
3587                {
[147167f]3588                        facet *fDeleteMarker=NULL;
[73809b]3589                        slAct = slHead;
3590                        if(slAct==NULL)
3591                        {
[63a2f8]3592                                printf("Zero length SLA\n");
[73809b]3593                                facet *fCopy;
[4199b3]3594                                fCopy = fAct;                           
[73809b]3595                                while(fCopy!=NULL)
3596                                {
[4199b3]3597                                        if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
[73809b]3598                                        {
3599                                                if(slAct==NULL)
3600                                                {
[4199b3]3601                                                        slAct = slAct->shallowCopy(*fCopy);//shallow copy constructor
3602                                                        slHead = slAct;                                                         
[73809b]3603                                                }
3604                                                else
3605                                                {
3606                                                        slAct->next = slAct->shallowCopy(*fCopy);
3607                                                        slAct = slAct->next;
3608                                                }
3609                                        }
3610                                        fCopy = fCopy->next;
3611                                }
[4199b3]3612                                break;  //WTF?
[73809b]3613                        }
3614                        /*Comparison starts here*/
3615                        while(slAct!=NULL)
3616                        {
3617                                removalOccured=FALSE;
[96e408]3618#ifndef NDEBUG
[11a7dc]3619        printf("Checking facet (");fAct->fNormal->show(1,1);printf(") against (");slAct->fNormal->show(1,1);printf(")\n");
[4199b3]3620#endif 
[73809b]3621                                if(areEqual2(fAct,slAct))
3622                                {
[147167f]3623                                        fDeleteMarker=slAct;
[73809b]3624                                        if(slAct==slHead)
[63a2f8]3625                                        {
[4199b3]3626//                                              fDeleteMarker=slHead;
3627//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
3628                                                slHead = slAct->next;                                           
3629//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
[73809b]3630                                                if(slHead!=NULL)
[63a2f8]3631                                                {
[73809b]3632                                                        slHead->prev = NULL;
[63a2f8]3633                                                }
3634                                                fDeleteMarker->shallowDelete();
[4199b3]3635                                                //delete fDeleteMarker;//NOTE this messes up fAct in noRevS!
3636//                                              printf("headUCN@enq=%i\n",slHead->getUCN());
[73809b]3637                                        }
3638                                        else if (slAct==slEnd)
3639                                        {
3640                                                slEnd=slEnd->prev;
3641                                                slEnd->next = NULL;
[63a2f8]3642                                                fDeleteMarker->shallowDelete();
[64b6c6]3643                                                delete(fDeleteMarker);
[4199b3]3644                                        }                                                               
[73809b]3645                                        else
3646                                        {
3647                                                slAct->prev->next = slAct->next;
3648                                                slAct->next->prev = slAct->prev;
[63a2f8]3649                                                fDeleteMarker->shallowDelete();
[64b6c6]3650                                                delete(fDeleteMarker);
[73809b]3651                                        }
3652                                        removalOccured=TRUE;
3653                                        gcone::lengthOfSearchList--;
[96e408]3654#ifndef NDEBUG
[63a2f8]3655printf("Removing (");fAct->fNormal->show(1,1);printf(") from list\n");
[73809b]3656#endif
[4199b3]3657                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
[73809b]3658                                }
[4199b3]3659                                slAct = slAct->next;                           
3660                        }//while(slAct!=NULL)
[73809b]3661                        if(removalOccured==FALSE)
3662                        {
3663                                facet *marker=slEnd;
3664                                slEnd->next = fAct->shallowCopy(*fAct);
3665                                slEnd = slEnd->next;
3666                                slEnd->prev=marker;
3667                                gcone::lengthOfSearchList++;
3668                        }
3669                        fAct = fAct->next;
[4199b3]3670//                      if(fDeleteMarker!=NULL)
3671//                      {
3672//                              fDeleteMarker->shallowDelete();
3673//                              delete(fDeleteMarker);
3674//                              fDeleteMarker=NULL;
3675//                      }
[73809b]3676                }
3677                else
3678                        fAct = fAct->next;
[4199b3]3679        }//while(fAct!=NULL)
3680       
[73809b]3681#ifdef gfanp
3682        gettimeofday(&end, 0);
3683        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
[4199b3]3684#endif 
3685//      printf("headUCN@enq=%i\n",slHead->getUCN());
[73809b]3686        return slHead;
3687}
3688
[4199b3]3689/**
[be7ab3f]3690* During flip2 every gc->baseRing gets two ringorder_a
3691* To avoid having an awful lot of those in the end we endow
3692* gc->baseRing by a suitable ring with (a,dp,C) and copy all
3693* necessary stuff over
3694* But first we will try to just do an inplace exchange and copying only the
3695* gc->gcBasis
3696*/
[ec6c52]3697void gcone::replaceDouble_ringorder_a_ByASingleOne()
[be7ab3f]3698{
3699        ring srcRing=currRing;
3700        ring replacementRing=rCopy0((ring)this->baseRing);
3701        /*We assume we have (a(),a(),dp) here*/
3702        omFree(replacementRing->order);
[590f813]3703        replacementRing->order =(int *)omAlloc0(4*sizeof(int/*64*/));
[be7ab3f]3704        omFree(replacementRing->block0);
[590f813]3705        replacementRing->block0=(int *)omAlloc0(4*sizeof(int/*64*/));
[be7ab3f]3706        omFree(replacementRing->block1);
[590f813]3707        replacementRing->block1=(int *)omAlloc0(4*sizeof(int/*64*/));
[be7ab3f]3708        omfree(replacementRing->wvhdl);
[590f813]3709        replacementRing->wvhdl =(int **)omAlloc0(4*sizeof(int/*64*/**));
[4199b3]3710                       
[590f813]3711        replacementRing->order[0]=ringorder_a/*64*/;
[be7ab3f]3712        replacementRing->block0[0]=1;
3713        replacementRing->block1[0]=replacementRing->N;
[4199b3]3714               
[be7ab3f]3715        replacementRing->order[1]=ringorder_dp;
3716        replacementRing->block0[1]=1;
3717        replacementRing->block1[1]=replacementRing->N;
[4199b3]3718       
[be7ab3f]3719        replacementRing->order[2]=ringorder_C;
[72e467]3720
[4199b3]3721        int64vec *ivw = this->getIntPoint(TRUE);//returns a reference
3722//      assert(this->ivIntPt); 
3723        int length=ivw->length();       
[590f813]3724        int/*64*/ *A=(int/*64*/ *)omAlloc0(length*sizeof(int/*64*/));
[be7ab3f]3725        for (int jj=0;jj<length;jj++)
[590f813]3726        {
[be7ab3f]3727                A[jj]=(*ivw)[jj];
[dd6b1d]3728                if((*ivw)[jj]>=INT_MAX) WarnS("A[jj] exceeds INT_MAX in gcone::replaceDouble_ringorder_a_ByASingleOne!\n");
[4199b3]3729        }       
3730        //delete ivw; //Not needed if this->getIntPoint(TRUE)
[be7ab3f]3731        replacementRing->wvhdl[0]=(int *)A;
3732        replacementRing->block1[0]=length;
3733        /*Finish*/
3734        rComplete(replacementRing);
3735        rChangeCurrRing(replacementRing);
3736        ideal temporaryGroebnerBasis=idrCopyR(this->gcBasis,this->baseRing);
3737        rDelete(this->baseRing);
3738        this->baseRing=rCopy(replacementRing);
[4199b3]3739        this->gcBasis=idCopy(temporaryGroebnerBasis);   
[be7ab3f]3740        /*And back to where we came from*/
[147167f]3741        rChangeCurrRing(srcRing);
3742        idDelete( (ideal*)&temporaryGroebnerBasis );
[4199b3]3743        rDelete(replacementRing);       
[be7ab3f]3744}
[72e467]3745
[1d9469]3746/** \brief Compute the gcd of two ints
3747 */
[f91fddc]3748static int64 int64gcd(const int64 &a, const int64 &b)
3749{
3750        int64 r, p=a, q=b;
3751        if(p < 0)
3752                p = -p;
3753        if(q < 0)
3754                q = -q;
3755        while(q != 0)
3756        {
3757                r = p % q;
3758                p = q;
3759                q = r;
3760        }
3761        return p;
3762}
[4199b3]3763       
3764static int intgcd(const int &a, const int &b)
3765{
3766        int r, p=a, q=b;
3767        if(p < 0)
3768                p = -p;
3769        if(q < 0)
3770                q = -q;
3771        while(q != 0)
3772        {
3773                r = p % q;
3774                p = q;
3775                q = r;
3776        }
3777        return p;
3778}
3779               
[1d9469]3780/** \brief Construct a dd_MatrixPtr from a cone's list of facets
[42660f]3781 * NO LONGER USED
[1d9469]3782 */
[4199b3]3783// inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
3784// {
3785//      facet *fAct;
3786//      fAct = gc.facetPtr;
3787//     
3788//      dd_MatrixPtr M;
3789//      dd_rowrange ddrows;
3790//      dd_colrange ddcols;
3791//      ddcols=(this->numVars)+1;
3792//      ddrows=this->numFacets;
3793//      dd_NumberType numb = dd_Integer;
3794//      M=dd_CreateMatrix(ddrows,ddcols);                       
3795//                     
3796//      int jj=0;
3797//     
3798//      while(fAct!=NULL)
3799//      {
3800//              int64vec *comp;
3801//              comp = fAct->getFacetNormal();
3802//              for(int ii=0;ii<this->numVars;ii++)
3803//              {
3804//                      dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
3805//              }
3806//              jj++;
3807//              delete comp;
3808//              fAct=fAct->next;                               
3809//      }                       
3810//      return M;
3811// }
3812               
[1d9469]3813/** \brief Write information about a cone into a file on disk
3814 *
3815 * This methods writes the information needed for the "second" method into a file.
3816 * The file's is divided in sections containing information on
3817 * 1) the ring
3818 * 2) the cone's Groebner Basis
3819 * 3) the cone's facets
3820 * Each line contains exactly one date
3821 * Each section starts with its name in CAPITALS
3822 */
[ec6c52]3823void gcone::writeConeToFile(const gcone &gc, bool usingIntPoints)
[1d9469]3824{
3825        int UCN=gc.UCN;
3826        stringstream ss;
3827        ss << UCN;
[4199b3]3828        string UCNstr = ss.str();               
3829                       
[17f35f6]3830        string prefix="/tmp/Singular/cone";
[4199b3]3831//      string prefix="./";     //crude hack -> run tests in separate directories
[bb503c7]3832        string suffix=".sg";
3833        string filename;
3834        filename.append(prefix);
3835        filename.append(UCNstr);
3836        filename.append(suffix);
[4199b3]3837               
3838//      int thisPID = getpid();
3839//      ss << thisPID;
3840//      string strPID = ss.str();
3841//      string prefix="./";
3842                                               
[bb503c7]3843        ofstream gcOutputFile(filename.c_str());
[73809b]3844        assert(gcOutputFile);
[bb503c7]3845        facet *fAct;
[4199b3]3846        fAct = gc.facetPtr;                     
[963eeb]3847        facet *f2Act;
3848        f2Act=fAct->codim2Ptr;
[4199b3]3849       
[be7ab3f]3850        char *ringString = rString(gc.baseRing);
[4199b3]3851                       
[1d9469]3852        if (!gcOutputFile)
3853        {
[11a7dc]3854                WerrorS("Error opening file for writing in writeConeToFile\n");
[1d9469]3855        }
3856        else
[4199b3]3857        {       
[963eeb]3858                gcOutputFile << "UCN" << endl;
3859                gcOutputFile << this->UCN << endl;
[4199b3]3860                gcOutputFile << "RING" << endl; 
[bec3b0]3861                gcOutputFile << ringString << endl;
3862                gcOutputFile << "GCBASISLENGTH" << endl;
[4199b3]3863                gcOutputFile << IDELEMS(gc.gcBasis) << endl;                   
3864                //Write this->gcBasis into file
3865                gcOutputFile << "GCBASIS" << endl;                             
[963eeb]3866                for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
[4199b3]3867                {                                       
[963eeb]3868                        gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
[4199b3]3869                }                               
3870                                       
3871                gcOutputFile << "FACETS" << endl;                                                               
3872               
[963eeb]3873                while(fAct!=NULL)
[4199b3]3874                {       
[f91fddc]3875                        const int64vec *iv=fAct->getRef2FacetNormal();
[4199b3]3876//                      iv=fAct->getRef2FacetNormal();//->getFacetNormal();
[963eeb]3877                        f2Act=fAct->codim2Ptr;
3878                        for (int ii=0;ii<iv->length();ii++)
[1d9469]3879                        {
[4199b3]3880//                              if (ii<iv->length()-1)
3881//                                      gcOutputFile << (*iv)[ii] << ",";
3882//                              else
3883//                                      gcOutputFile << (*iv)[ii] << " ";
[64b6c6]3884                                gcOutputFile << (*iv)[ii];
3885                                (ii<iv->length()-1) ? gcOutputFile << "," : gcOutputFile << " ";
[1d9469]3886                        }
[4199b3]3887                        //delete iv;
[963eeb]3888                        while(f2Act!=NULL)
[1d9469]3889                        {
[590f813]3890                                const int64vec *iv2;
[64b6c6]3891                                iv2=f2Act->getRef2FacetNormal();
[963eeb]3892                                for(int jj=0;jj<iv2->length();jj++)
3893                                {
[4199b3]3894//                                      if (jj<iv2->length()-1)
3895//                                              gcOutputFile << (*iv2)[jj] << ",";
3896//                                      else
3897//                                              gcOutputFile << (*iv2)[jj] << " ";
[64b6c6]3898                                        gcOutputFile << (*iv2)[jj];
3899                                        (jj<iv2->length()-1) ? gcOutputFile << "," : gcOutputFile << " ";
[963eeb]3900                                }
3901                                f2Act = f2Act->next;
[1d9469]3902                        }
[963eeb]3903                        gcOutputFile << endl;
3904                        fAct=fAct->next;
[4199b3]3905                }                       
[1d9469]3906        gcOutputFile.close();
3907        }
[f91fddc]3908        delete [] ringString;
[4199b3]3909                       
3910}//writeConeToFile(gcone const &gc)
3911               
[1d9469]3912/** \brief Reads a cone from a file identified by its number
[be7ab3f]3913* ||depending on whether flip or flip2 is used, switch the flag flipFlag
3914* ||defaults to 0 => flip
3915* ||1 => flip2
3916*/
[ec6c52]3917void gcone::readConeFromFile(int UCN, gcone *gc)
[1d9469]3918{
[4199b3]3919        //int UCN=gc.UCN;
[963eeb]3920        stringstream ss;
3921        ss << UCN;
[bec3b0]3922        string UCNstr = ss.str();
3923        int gcBasisLength=0;
[4199b3]3924        size_t found;   //used for find_first_(not)_of
[bb503c7]3925
[17f35f6]3926        string prefix="/tmp/Singular/cone";
[e98068]3927        string suffix=".sg";
3928        string filename;
3929        filename.append(prefix);
3930        filename.append(UCNstr);
3931        filename.append(suffix);
[4199b3]3932                                       
[e98068]3933        ifstream gcInputFile(filename.c_str(), ifstream::in);
[4199b3]3934       
3935        ring saveRing=currRing; 
3936        //Comment the following if you uncomment the if(line=="RING") part below
3937//      rChangeCurrRing(gc->baseRing);
3938       
[bec3b0]3939        while( !gcInputFile.eof() )
3940        {
[42660f]3941                string line;
[bec3b0]3942                getline(gcInputFile,line);
[a23a297]3943                if(line=="RING")
[2e06300]3944                {
[1f3450]3945                        getline(gcInputFile,line);
3946                        found = line.find("a(");
3947                        line.erase(0,found+2);
3948                        string strweight;
3949                        strweight=line.substr(0,line.find_first_of(")"));
[4199b3]3950                                               
3951                        int64vec *iv=new int64vec(this->numVars);//                     
[1f3450]3952                        for(int ii=0;ii<this->numVars;ii++)
3953                        {
3954                                string weight;
[dd6b1d]3955                                weight=line.substr(0,line.find_first_of(",)"));
3956                                char *w=new char[weight.size()+1];
3957                                strcpy(w,weight.c_str());
[4199b3]3958                                (*iv)[ii]=atol(w/*weight.c_str()*/);//Better to long. Weight bound in Singular:2147483647
[dd6b1d]3959                                delete[] w;
[1f3450]3960                                line.erase(0,line.find_first_of(",)")+1);
3961                        }
[be7ab3f]3962                        found = line.find("a(");
[4199b3]3963                       
[1f3450]3964                        ring newRing;
[590f813]3965                        if(currRing->order[0]!=ringorder_a/*64*/)
[1f3450]3966                        {
[be7ab3f]3967                                        newRing=rCopyAndAddWeight(currRing,iv);
[1f3450]3968                        }
3969                        else
[4199b3]3970                        {       
[be7ab3f]3971                                        newRing=rCopy0(currRing);
3972                                        int length=this->numVars;
3973                                        int *A=(int *)omAlloc0(length*sizeof(int));
3974                                        for(int jj=0;jj<length;jj++)
3975                                        {
3976                                                A[jj]=(*iv)[jj];
3977                                        }
3978                                        omFree(newRing->wvhdl[0]);
3979                                        newRing->wvhdl[0]=(int*)A;
3980                                        newRing->block1[0]=length;
[1f3450]3981                        }
3982                        delete iv;
3983                        rComplete(newRing);
3984                        gc->baseRing=rCopy(newRing);
3985                        rDelete(newRing);
3986                        rComplete(gc->baseRing);
3987                        if(currRing!=gc->baseRing)
3988                                rChangeCurrRing(gc->baseRing);
[a23a297]3989                }
[4199b3]3990               
[bec3b0]3991                if(line=="GCBASISLENGTH")
3992                {
[42660f]3993                        string strGcBasisLength;
[bec3b0]3994                        getline(gcInputFile, line);
[e98068]3995                        strGcBasisLength = line;
[dd6b1d]3996                        char *s=new char[strGcBasisLength.size()+1];
3997                        strcpy(s,strGcBasisLength.c_str());
3998                        int size=atoi(s/*strGcBasisLength.c_str()*/);
3999                        delete[] s;
[4199b3]4000                        gcBasisLength=size;             
[6123d46]4001                        gc->gcBasis=idInit(size,1);
[bec3b0]4002                }
4003                if(line=="GCBASIS")
4004                {
4005                        for(int jj=0;jj<gcBasisLength;jj++)
4006                        {
4007                                getline(gcInputFile,line);
[4199b3]4008                                //magically convert strings into polynomials
4009                                //polys.cc:p_Read
4010                                //check until first occurance of + or -
4011                                //data or c_str
4012//                              poly strPoly;//=pInit();//Ought to be inside the while loop, but that will eat your memory
4013                                poly resPoly=pInit();   //The poly to be read in               
[e98068]4014                                while(!line.empty())
[bec3b0]4015                                {
[4199b3]4016                                        poly strPoly;//=pInit();
4017                                       
[42660f]4018                                        string strMonom, strCoeff, strCoeffNom, strCoeffDenom;
[4199b3]4019                                        bool hasCoeffInQ = FALSE;       //does the polynomial have rational coeff?
4020                                        bool hasNegCoeff = FALSE;       //or a negative one?
4021                                        found = line.find_first_of("+-");       //get the first monomial
[e98068]4022                                        string tmp;
4023                                        tmp=line[found];
[4199b3]4024//                                      if(found!=0 && (tmp.compare("-")==0) )
4025//                                              hasNegCoeff = TRUE;     //We have a coeff < 0
[e98068]4026                                        if(found==0)
4027                                        {
4028                                                if(tmp.compare("-")==0)
4029                                                        hasNegCoeff = TRUE;
[4199b3]4030                                                line.erase(0,1);        //remove leading + or -
4031                                                found = line.find_first_of("+-");       //adjust found
[e98068]4032                                        }
4033                                        strMonom = line.substr(0,found);
4034                                        line.erase(0,found);
[dd6b1d]4035                                        number nCoeff=nInit(1);
4036                                        number nCoeffNom=nInit(1);
4037                                        number nCoeffDenom=nInit(1);
[e98068]4038                                        found = strMonom.find_first_of("/");
[4199b3]4039                                        if(found!=string::npos) //i.e. "/" exists in strMonom
[e98068]4040                                        {
4041                                                hasCoeffInQ = TRUE;
[4199b3]4042                                                strCoeffNom=strMonom.substr(0,found);                                           
[f07b38c]4043                                                strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890",found+1));
[d2cd515]4044                                                strMonom.erase(0,found);
[4199b3]4045                                                strMonom.erase(0,strMonom.find_first_not_of("1234567890/"));   
[dd6b1d]4046                                                char *Nom=new char[strCoeffNom.size()+1];
4047                                                char *Denom=new char[strCoeffDenom.size()+1];
4048                                                strcpy(Nom,strCoeffNom.c_str());
[4199b3]4049                                                strcpy(Denom,strCoeffDenom.c_str());           
[dd6b1d]4050                                                nRead(Nom/*strCoeffNom.c_str()*/, &nCoeffNom);
4051                                                nRead(Denom/*strCoeffDenom.c_str()*/, &nCoeffDenom);
4052                                                delete[] Nom;
4053                                                delete[] Denom;
[e98068]4054                                        }
4055                                        else
[bec3b0]4056                                        {
[e98068]4057                                                found = strMonom.find_first_not_of("1234567890");
[4199b3]4058                                                strCoeff = strMonom.substr(0,found);                                           
[e98068]4059                                                if(!strCoeff.empty())
4060                                                {
[dd6b1d]4061                                                        char *coeff = new char[strCoeff.size()+1];
4062                                                        strcpy(coeff, strCoeff.c_str());
4063                                                        nRead(coeff/*strCoeff.c_str()*/,&nCoeff);
4064                                                        delete[] coeff;
[e98068]4065                                                }
4066                                        }
[dd6b1d]4067                                        char* monom = new char[strMonom.size()+1];
[4199b3]4068                                        strcpy(monom, strMonom.c_str());                                               
4069                                        p_Read(monom,strPoly,currRing); //strPoly:=monom
4070                                        delete[] monom; 
[e98068]4071                                        switch (hasCoeffInQ)
4072                                        {
4073                                                case TRUE:
4074                                                        if(hasNegCoeff)
[d2cd515]4075                                                                nCoeffNom=nNeg(nCoeffNom);
[e98068]4076                                                        pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom));
[d2cd515]4077                                                        break;
[e98068]4078                                                case FALSE:
4079                                                        if(hasNegCoeff)
[4199b3]4080                                                                nCoeff=nNeg(nCoeff);                                                   
[e98068]4081                                                        if(!nIsOne(nCoeff))
[bec3b0]4082                                                        {
[e98068]4083                                                                pSetCoeff(strPoly, nCoeff );
[bec3b0]4084                                                        }
[d2cd515]4085                                                        break;
[bec3b0]4086                                        }
[4199b3]4087                                                //pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1???
[e98068]4088                                        if(pIsConstantComp(resPoly))
[dd6b1d]4089                                        {
4090                                                resPoly=pCopy(strPoly);
4091                                                pDelete(&strPoly);
4092                                        }
[e98068]4093                                        else
[dd6b1d]4094                                        {
[4199b3]4095//                                              poly tmp=pAdd(pCopy(resPoly),strPoly);//foo is destroyed
4096//                                              pDelete(&resPoly);
4097//                                              resPoly=tmp;
4098//                                              pDelete(&tmp);
4099                                                resPoly=pAdd(resPoly,strPoly);//pAdd = p_Add_q, destroys args
[dd6b1d]4100                                        }
[4199b3]4101                                        /*if(nCoeff!=NULL)                                     
4102                                                nDelete(&nCoeff);*/ //NOTE This may cause a crash on certain examples...
[42660f]4103                                        nDelete(&nCoeffNom);
4104                                        nDelete(&nCoeffDenom);
[4199b3]4105                                }//while(!line.empty())                 
[f07b38c]4106                                gc->gcBasis->m[jj]=pCopy(resPoly);
[4199b3]4107                                pDelete(&resPoly);      //reset
[bec3b0]4108                        }
[4199b3]4109//                      break;
4110                }//if(line=="GCBASIS") 
[706b89]4111                if(line=="FACETS")
4112                {
4113                        facet *fAct=gc->facetPtr;
[63a2f8]4114                        while(fAct!=NULL)
[706b89]4115                        {
4116                                getline(gcInputFile,line);
4117                                found = line.find("\t");
4118                                string normalString=line.substr(0,found);
[590f813]4119                                int64vec *fN = new int64vec(this->numVars);
[706b89]4120                                for(int ii=0;ii<this->numVars;ii++)
4121                                {
4122                                        string component;
4123                                        found = normalString.find(",");
4124                                        component=normalString.substr(0,found);
[dd6b1d]4125                                        char *sComp = new char[component.size()+1];
4126                                        strcpy(sComp,component.c_str());
4127                                        (*fN)[ii]=atol(sComp/*component.c_str()*/);
4128                                        delete[] sComp;
[706b89]4129                                        normalString.erase(0,found+1);
4130                                }
4131                                /*Only the following line needs to be commented out if you decide not to delete fNormals*/
[4199b3]4132//                              fAct->setFacetNormal(fN);
[706b89]4133                                delete(fN);
[4199b3]4134                                fAct = fAct->next;      //Booh, this is ugly
4135                        }                       
4136                        break; //NOTE Must always be in the last if-block!
[706b89]4137                }
[4199b3]4138        }//while(!gcInputFile.eof())   
[963eeb]4139        gcInputFile.close();
[bec3b0]4140        rChangeCurrRing(saveRing);
[7c0ec6]4141}
[4199b3]4142       
[73809b]4143
[5ece85]4144/** \brief Sort the rays of a facet lexicographically
4145*/
[4199b3]4146// void gcone::sortRays(gcone *gc)
4147// {
4148//      facet *fAct;
4149//      fAct = this->facetPtr->codim2Ptr;
4150//      while(fAct->next!=NULL)
4151//      {
4152//              if(fAct->fNormal->compare(fAct->fNormal->next)==-1
4153//      }
4154// }
[5ece85]4155
[0d278cd]4156/** \brief Gather the output
4157* List of lists
[1f3450]4158* If heuristic==1 readConeFromFile() is called once more on every cone. This may slow down the computation but it also
4159* allows us to rDelete(gcDel->baseRing) and the such in gcone::noRevS.
[b45f583]4160*\param *gc Pointer to gcone, preferably gcRoot ;-)
[276932]4161*\param n the number of cones as determined by gcRoot->getCounter()
[b45f583]4162*
[0d278cd]4163*/
[73809b]4164lists lprepareResult(gcone *gc, const int n)
[005d00a]4165{
4166        gcone *gcAct;
[4199b3]4167        gcAct = gc;     
[005d00a]4168        facet *fAct;
[0d278cd]4169        fAct = gc->facetPtr;
[4199b3]4170       
[0d278cd]4171        lists res=(lists)omAllocBin(slists_bin);
[4199b3]4172        res->Init(n);   //initialize to store n cones
[0d278cd]4173        for(int ii=0;ii<n;ii++)
[005d00a]4174        {
[4199b3]4175                if(gfanHeuristic==1)// && gcAct->getUCN()>1)
[276932]4176                {
[4199b3]4177                        gcAct->readConeFromFile(gcAct->getUCN(),gcAct); 
4178//                      rChangeCurrRing(gcAct->getBaseRing());//NOTE memleak?
[147167f]4179                }
[4199b3]4180                rChangeCurrRing(gcAct->getRef2BaseRing());     
[0d278cd]4181                res->m[ii].rtyp=LIST_CMD;
4182                lists l=(lists)omAllocBin(slists_bin);
[e2e3ad]4183                l->Init(6);
[0d278cd]4184                l->m[0].rtyp=INT_CMD;
4185                l->m[0].data=(void*)gcAct->getUCN();
[b45f583]4186                l->m[1].rtyp=IDEAL_CMD;
[f07b38c]4187                /*The following is necessary for leaves in the tree of cones
[4199b3]4188                * Since we don't use them in the computation and gcBasis is
[f07b38c]4189                * set to (poly)NULL in noRevS we need to get this back here.
4190                */
[4199b3]4191//              if(gcAct->gcBasis->m[0]==(poly)NULL)
4192//              if(gfanHeuristic==1 && gcAct->getUCN()>1)
4193//                      gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
4194//              ring saveRing=currRing;
4195//              ring tmpRing=gcAct->getBaseRing;
4196//              rChangeCurrRing(tmpRing);
4197//              l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
4198//              l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());//NOTE memleak?
[147167f]4199                l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getRef2BaseRing());
[4199b3]4200//              rChangeCurrRing(saveRing);
[f07b38c]4201
[b45f583]4202                l->m[2].rtyp=INTVEC_CMD;
[4199b3]4203                int64vec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));//NOTE memleak?
[590f813]4204                l->m[2].data=(void*)iv64Copy(&iv);
[4199b3]4205               
[b45f583]4206                l->m[3].rtyp=LIST_CMD;
4207                        lists lCodim2List = (lists)omAllocBin(slists_bin);
[4199b3]4208                        lCodim2List->Init(gcAct->numFacets); 
4209                        fAct = gcAct->facetPtr;//fAct->codim2Ptr;
[b45f583]4210                        int jj=0;
4211                        while(fAct!=NULL && jj<gcAct->numFacets)
4212                        {
4213                                lCodim2List->m[jj].rtyp=INTVEC_CMD;
[590f813]4214                                int64vec ivC2=(gcAct->f2M(gcAct,fAct,2));
4215                                lCodim2List->m[jj].data=(void*)iv64Copy(&ivC2);
[b45f583]4216                                jj++;
4217                                fAct = fAct->next;
[4199b3]4218                        }               
4219                l->m[3].data=(void*)lCodim2List;               
[231da1]4220                l->m[4].rtyp=INTVEC_CMD/*RING_CMD*/;
[4199b3]4221                l->m[4].data=(void*)(gcAct->getIntPoint/*BaseRing*/());                 
[e2e3ad]4222                l->m[5].rtyp=INT_CMD;
4223                l->m[5].data=(void*)gcAct->getPredUCN();
[0d278cd]4224                res->m[ii].data=(void*)l;
[005d00a]4225                gcAct = gcAct->next;
[4199b3]4226        }       
[0d278cd]4227        return res;
[005d00a]4228}
[4199b3]4229
4230/** Convert gc->gcRays into an intvec in order to be used with bbcone stuff*/
4231intvec *gcRays2Intmat(gcone *gc)
4232{
4233  int r = gc->numRays;
4234  int c = gc->numVars;  //Spalten
4235  intvec *res = new intvec(r,c,(int)0);
4236 
4237  int offset=0;
4238  for(int ii=0;ii<gc->numRays;ii++)
4239  {
4240    int64vec *ivTmp=iv64Copy(gc->gcRays[ii]);
[1f637e]4241    for(int jj=0;jj<(currRing->N);jj++)
[4199b3]4242      (*res)[offset+jj]=(int)(*ivTmp)[jj];
[1f637e]4243    offset += (currRing->N);
[4199b3]4244    delete ivTmp;
4245  }
4246  return res;
4247}
4248
4249/** \brief Put stuff in gfanlib's datatype gfan::ZFan
4250*/
4251void prepareGfanLib(gcone *gc, gfan::ZFan *fan)
4252{
4253  using namespace gfan;
4254  int ambientDimension = gc->numVars;
4255  gcone *gcAct;
4256  gcAct = gc;
4257
4258  //Iterate over all cones and adjoin to PolyhedralFan
4259   while(gcAct!=NULL)
4260  {     
4261    intvec *rays=gcRays2Intmat(gcAct);
4262    ZMatrix zm = intmat2ZMatrix(rays);
4263    delete rays;
4264    ZCone *zc = new ZCone();
4265    *zc = ZCone::givenByRays(zm, gfan::ZMatrix(0, zm.getWidth()));
4266//     delete &zm;
[d503ac]4267    zc->canonicalize();//As per Anders' hint
[4199b3]4268    fan->insert(*zc);
4269//     delete zc;
4270    gcAct=gcAct->next;
4271  }
4272}
4273
[005d00a]4274/** \brief Write facets of a cone into a matrix
[4199b3]4275* Takes a pointer to a facet as 2nd arg
[b45f583]4276* f should always point to gc->facetPtr
4277* param n is used to determine whether it operates in codim 1 or 2
[590f813]4278* We have to cast the int64vecs to int64vec due to issues with list structure
[005d00a]4279*/
[590f813]4280inline int64vec gcone::f2M(gcone *gc, facet *f, int n)
[005d00a]4281{
4282        facet *fAct;
[4199b3]4283        int64vec *res;//=new int64vec(this->numVars);   
4284//      int codim=n;
4285//      int bound;
4286//      if(f==gc->facetPtr)
[b45f583]4287        if(n==1)
4288        {
[590f813]4289                int64vec *m1Res=new int64vec(gc->numFacets,gc->numVars,0);
4290                res = iv64Copy(m1Res);
[005d00a]4291                fAct = gc->facetPtr;
[dcf8b4b]4292                delete m1Res;
[4199b3]4293//              bound = gc->numFacets*(this->numVars);         
[b45f583]4294        }
[005d00a]4295        else
4296        {
[b45f583]4297                fAct = f->codim2Ptr;
[590f813]4298                int64vec *m2Res = new int64vec(f->numCodim2Facets,gc->numVars,0);
4299                res = iv64Copy(m2Res);
[4199b3]4300                delete m2Res;   
4301//              bound = fAct->numCodim2Facets*(this->numVars);
[b45f583]4302
[005d00a]4303        }
[0d278cd]4304        int ii=0;
[4199b3]4305        while(fAct!=NULL )//&& ii < bound )
[005d00a]4306        {
[590f813]4307                const int64vec *fNormal;
[4199b3]4308                fNormal = fAct->getRef2FacetNormal();//->getFacetNormal();
[b45f583]4309                for(int jj=0;jj<this->numVars;jj++)
4310                {
[4199b3]4311                        (*res)[ii]=(int)(*fNormal)[jj];//This is ugly and prone to overflow
[005d00a]4312                        ii++;
4313                }
4314                fAct = fAct->next;
[4199b3]4315        }       
[b45f583]4316        return *res;
[005d00a]4317}
[2c6535]4318
[b7510d]4319int gcone::counter=0;
[d2cd515]4320int gfanHeuristic;
[72e467]4321int gcone::lengthOfSearchList;
4322int gcone::maxSize;
4323dd_MatrixPtr gcone::dd_LinealitySpace;
[590f813]4324int64vec *gcone::hilbertFunction;
[42660f]4325#ifdef gfanp
[4199b3]4326// int gcone::lengthOfSearchList=0;
[42660f]4327float gcone::time_getConeNormals;
4328float gcone::time_getCodim2Normals;
[4c4a80]4329float gcone::t_getExtremalRays;
[73809b]4330float gcone::t_ddPolyh;
[42660f]4331float gcone::time_flip;
[be7ab3f]4332float gcone::time_flip2;
[85583b]4333float gcone::t_areEqual;
[2e06300]4334float gcone::t_markings;
4335float gcone::t_dd;
[11a7dc]4336float gcone::t_kStd=0;
[42660f]4337float gcone::time_enqueue;
4338float gcone::time_computeInv;
4339float gcone::t_ddMC;
4340float gcone::t_mI;
4341float gcone::t_iP;
[73809b]4342float gcone::t_isParallel;
[4c4a80]4343unsigned gcone::parallelButNotEqual=0;
4344unsigned gcone::numberOfFacetChecks=0;
[42660f]4345#endif
[5151ac]4346int gcone::numVars;
[72e467]4347bool gcone::hasHomInput=FALSE;
[590f813]4348int64vec *gcone::ivZeroVector;
[4199b3]4349// ideal gfan(ideal inputIdeal, int h)
4350/** Main routine
[305057]4351 * The first and second parameter are mandatory. The third (and maybe fourth) parameter is for Janko :)
4352 */
[4199b3]4353#ifndef USE_ZFAN
4354lists grfan(ideal inputIdeal, int h, bool singleCone=FALSE)
4355#else
4356gfan::ZFan* grfan(ideal inputIdeal, int h, bool singleCone=FALSE)
4357#endif
[2c6535]4358{
[4199b3]4359        lists lResList; //this is the object we return 
[1f637e]4360        gfan::ZFan *zResFan = new gfan::ZFan((currRing->N));
[4199b3]4361       
[197a82]4362        if(rHasGlobalOrdering(currRing))
[4199b3]4363        {       
[1f637e]4364//              int numvar = (currRing->N);
[4c4a80]4365                gfanHeuristic = h;
[4199b3]4366               
[4c4a80]4367                enum searchMethod {
4368                        reverseSearch,
4369                        noRevS
4370                };
[4199b3]4371       
[4c4a80]4372                searchMethod method;
4373                method = noRevS;
[4199b3]4374       
4375                ring inputRing=currRing;        // The ring the user entered
4376//              ring rootRing;                  // The ring associated to the target ordering
[276932]4377
[4c4a80]4378                dd_set_global_constants();
4379                if(method==noRevS)
[913422]4380                {
[4c4a80]4381                        gcone *gcRoot = new gcone(currRing,inputIdeal);
4382                        gcone *gcAct;
4383                        gcAct = gcRoot;
[1f637e]4384                        gcone::numVars=(currRing->N);
4385                        //gcAct->numVars=(currRing->N);//NOTE is now static
[4c4a80]4386                        gcAct->getGB(inputIdeal);
4387                        /*Check whether input is homogeneous
4388                        if TRUE each facet intersects the positive orthant, so we don't need the
4389                        flippability test in getConeNormals & getExtremalRays
4390                        */
[4199b3]4391                        if(idHomIdeal(gcAct->gcBasis,NULL))//disabled for tests
[2f72c55]4392                        {
[4c4a80]4393                                gcone::hasHomInput=TRUE;
[4199b3]4394//                              gcone::hilbertFunction=hHstdSeries(inputIdeal,NULL,NULL,NULL,currRing);
[2f72c55]4395                        }
[b3bb13]4396                        else
4397                        {
[1f637e]4398                                gcone::ivZeroVector = new int64vec((currRing->N));
4399                                for(int ii=0;ii<(currRing->N);ii++)
[b3bb13]4400                                        (*gcone::ivZeroVector)[ii]=0;
4401                        }
[4199b3]4402
[ec6c52]4403                        if(isMonomial(gcAct->gcBasis))
[4199b3]4404                        {//FIXME
[4c4a80]4405                                WerrorS("Monomial input - terminating");
[b3bb13]4406                                dd_free_global_constants();
[4199b3]4407                                //This is filthy
4408                                goto pointOfNoReturn;                           
4409                        }                       
[dcf8b4b]4410                        gcAct->getConeNormals(gcAct->gcBasis);
[be7ab3f]4411                        gcone::dd_LinealitySpace = gcAct->computeLinealitySpace();
4412                        gcAct->getExtremalRays(*gcAct);
[4199b3]4413                        if(singleCone==FALSE)//Is Janko here?
4414                        {//Compute the whole fan
4415                          gcAct->noRevS(*gcAct);        //Here we go!
4416                        }                           
4417                        //Switch back to the ring the computation was started in
[ae7cad0]4418                        rChangeCurrRing(inputRing);
[4199b3]4419                        //res=gcAct->gcBasis;
4420                        //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
4421#ifndef USE_ZFAN
[4c4a80]4422                        lResList=lprepareResult(gcRoot,gcRoot->getCounter());
[4199b3]4423#else
4424                        prepareGfanLib(gcRoot,zResFan);
4425#endif
[4c4a80]4426                        /*Cleanup*/
[4199b3]4427                        gcone *gcDel;   
[4c4a80]4428                        gcDel = gcRoot;
4429                        gcAct = gcRoot;
4430                        while(gcAct!=NULL)
4431                        {
4432                                gcDel = gcAct;
4433                                gcAct = gcAct->next;
[4199b3]4434//                              delete gcDel;
[4c4a80]4435                        }
[4199b3]4436                }//method==noRevS
[4c4a80]4437                dd_FreeMatrix(gcone::dd_LinealitySpace);
4438                dd_free_global_constants();
[4199b3]4439        }//rHasGlobalOrdering
[197a82]4440        else
4441        {
[4199b3]4442                //Simply return an empty list
[11a7dc]4443                WerrorS("Ring has non-global ordering.\nThis function requires your current ring to be endowed with a global ordering.\n Now terminating!");
[4199b3]4444//              gcone *gcRoot=new gcone();
4445//              gcone *gcPtr = gcRoot;
4446//              for(int ii=0;ii<10000;ii++)
4447//              {
4448//                      gcPtr->setBaseRing(currRing);
4449//                      facet *fPtr=gcPtr->facetPtr=new facet();
4450//                      for(int jj=0;jj<5;jj++)
4451//                      {
[1f637e]4452//                              int64vec *iv=new int64vec((currRing->N));
[4199b3]4453//                              fPtr->setFacetNormal(iv);                               
4454//                              delete(iv);
4455//                              fPtr->next=new facet();
4456//                              fPtr=fPtr->next;
4457//                      }
4458//                      gcPtr->next=new gcone();
4459//                      gcPtr->next->prev=gcPtr;
4460//                      gcPtr=gcPtr->next;                     
4461//              }
4462//              gcPtr=gcRoot;
4463//              while(gcPtr!=NULL)
4464//              {
4465//                      gcPtr=gcPtr->next;
4466// //                   delete(gcPtr->prev);
4467//              }
[11a7dc]4468                goto pointOfNoReturn;
[197a82]4469        }
[f07b38c]4470        /*Return result*/
[42660f]4471#ifdef gfanp
[1f3450]4472        cout << endl << "t_getConeNormals:" << gcone::time_getConeNormals << endl;
[64b6c6]4473        /*cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl;
4474        cout << "  t_ddMC:" << gcone::t_ddMC << endl;
4475        cout << "  t_mI:" << gcone::t_mI << endl;
4476        cout << "  t_iP:" << gcone::t_iP << endl;*/
[4c4a80]4477        cout << "t_getExtremalRays:" << gcone::t_getExtremalRays << endl;
[73809b]4478        cout << "  t_ddPolyh:" << gcone::t_ddPolyh << endl;
[42660f]4479        cout << "t_Flip:" << gcone::time_flip << endl;
[2e06300]4480        cout << "  t_markings:" << gcone::t_markings << endl;
4481        cout << "  t_dd:" << gcone::t_dd << endl;
[85583b]4482        cout << "  t_kStd:" << gcone::t_kStd << endl;
[be7ab3f]4483        cout << "t_Flip2:" << gcone::time_flip2 << endl;
4484        cout << "  t_dd:" << gcone::t_dd << endl;
4485        cout << "  t_kStd:" << gcone::t_kStd << endl;
[42660f]4486        cout << "t_computeInv:" << gcone::time_computeInv << endl;
4487        cout << "t_enqueue:" << gcone::time_enqueue << endl;
[85583b]4488        cout << "  t_areEqual:" <<gcone::t_areEqual << endl;
[73809b]4489        cout << "t_isParallel:" <<gcone::t_isParallel << endl;
[4c4a80]4490        cout << endl;
4491        cout << "Checked " << gcone::numberOfFacetChecks << " Facets" << endl;
4492        cout << " out of which there were " << gcone::parallelButNotEqual << " parallel but not equal." << endl;
[42660f]4493#endif
[11a7dc]4494        printf("Maximum lenght of list of facets: %i", gcone::maxSize);
4495pointOfNoReturn:
[4199b3]4496#ifndef USE_ZFAN
[005d00a]4497        return lResList;
[4199b3]4498#else
4499        return zResFan;
4500#endif
[d3398c]4501}
[197a82]4502
[305057]4503/** Compute a single Gröbner cone by specifying an ideal and a weight vector.
4504 * NOTE: We do NOT check whether the vector is from the relative interior of the cone.
4505 * That is upon the user to assure.
4506 */
[4199b3]4507// lists grcone_by_intvec(ideal inputIdeal)
4508// {
4509//   if( (rRingOrder_t)currRing->order[0] == ringorder_wp)
4510//   {
4511//     lists lResList;
4512//     lResList=grfan(inputIdeal, 0, TRUE);
4513//   }
4514//   else
4515//     WerrorS("Need wp ordering");
4516// }
[34fcf8]4517#endif
Note: See TracBrowser for help on using the repository browser.