source: git/kernel/gfan.cc @ 2e06300

spielwiese
Last change on this file since 2e06300 was 2e06300, checked in by Martin Monerjan, 14 years ago
Various attempts on preprocessing gcone::preprocessInequalities is the one to proceed with interiorPoint2 as a crude attempt on computing the said Possibility to remove strictly positive rows in getConeNormals - no speedup worth mentioning Several minor tweaks - note that getCodim2Normals actually seems to compute the extremal rays. However that does not do any harm since comparison works that way as well. Might only come cheaper to get cddlib to return the actual normals. But who knows with cddlib... Timestamps with gettimeofday for crude profiling dd_FindRelativeInterior -> dd_LPSolve Ring construction in readConeFromFile git-svn-id: file:///usr/local/Singular/svn/trunk@12395 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 81.2 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
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 $
6$Id$
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_GFAN
12
13#include "kstd1.h"
14#include "kutil.h"      //ksCreateSpoly
15// #include "intvec.h"
16#include "polys.h"
17#include "ideals.h"
18#include "kmatrix.h"
19#include "GMPrat.h"
20//#include "fast_maps.h"        //Mapping of ideals
21// #include "maps.h"
22// #include "ring.h"    //apparently not needed
23// #include "structs.h"
24#include "../Singular/lists.h"
25#include "prCopy.h"     //Needed for idrCopyR
26#include <iostream>
27// #include <bitset>
28#include <fstream>      //read-write cones to files
29// #include <gmp.h>
30#include <string>
31#include <sstream>
32#include <time.h>
33#include <sys/time.h>
34//#include <gmpxx.h>
35
36/*DO NOT REMOVE THIS*/
37#ifndef GMPRATIONAL
38#define GMPRATIONAL
39#endif
40
41//Hacks for different working places
42#define p800
43
44#ifdef UNI
45#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
46#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
47#endif
48
49#ifdef HOME
50#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
51#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
52#endif
53
54#ifdef ITWM
55#include "/u/slg/monerjan/cddlib/include/setoper.h"
56#include "/u/slg/monerjan/cddlib/include/cdd.h"
57#include "/u/slg/monerjan/cddlib/include/cddmp.h"
58#endif
59
60#ifdef p800
61#include "../../cddlib/include/setoper.h"
62#include "../../cddlib/include/cdd.h"
63#include "../../cddlib/include/cddmp.h"
64#endif
65
66#ifndef gfan_DEBUG
67// #define gfan_DEBUG
68#ifndef gfan_DEBUGLEVEL
69#define gfan_DEBUGLEVEL 1
70#endif
71#endif
72#define gfanp
73#include <gfan.h>
74using namespace std;
75
76/**
77*\brief Class facet
78*       Implements the facet structure as a linked list
79*
80*/
81               
82/** \brief The default constructor for facets
83*/
84facet::facet()                 
85{
86        // Pointer to next facet.  */
87        /* Defaults to NULL. This way there is no need to check explicitly */
88        this->fNormal=NULL;
89        this->interiorPoint=NULL;               
90        this->UCN=0;
91        this->codim2Ptr=NULL;
92        this->codim=1;          //default for (codim-1)-facets
93        this->numCodim2Facets=0;
94        this->flipGB=NULL;
95        this->isIncoming=FALSE;
96        this->next=NULL;
97        this->prev=NULL;               
98        this->flipRing=NULL;    //the ring on the other side
99        this->isFlippable=FALSE;
100}
101               
102/** \brief Constructor for facets of codim >= 2
103* Note that as of now the code of the constructors is only for facets and codim2-faces. One
104* could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
105*/
106facet::facet(int const &n)
107{
108        this->fNormal=NULL;
109        this->interiorPoint=NULL;                       
110        this->UCN=0;
111        this->codim2Ptr=NULL;
112        if(n>1)
113        {
114                this->codim=n;
115        }//NOTE Handle exception here!                 
116        this->numCodim2Facets=0;
117        this->flipGB=NULL;
118        this->isIncoming=FALSE;
119        this->next=NULL;
120        this->prev=NULL;
121        this->flipRing=NULL;
122        this->isFlippable=FALSE;
123}
124               
125/** \brief The copy constructor
126* By default only copies the fNormal, f2Normals and UCN
127*/
128facet::facet(const facet& f)
129{
130        this->fNormal=ivCopy(f.fNormal);
131        this->UCN=f.UCN;
132        this->isFlippable=f.isFlippable;
133        facet* f2Copy;
134        f2Copy=f.codim2Ptr;
135        facet* f2Act;
136        f2Act=this->codim2Ptr;
137        while(f2Copy!=NULL)
138        {
139                if(f2Act==NULL)
140                {
141                        f2Act=new facet(2);
142                        this->codim2Ptr=f2Act;                 
143                }
144                else
145                {
146                        facet* marker;
147                        marker = f2Act;
148                        f2Act->next = new facet(2);
149                        f2Act = f2Act->next;
150                        f2Act->prev = marker;
151                }
152                intvec *f2Normal;
153                f2Normal = f2Copy->getFacetNormal();
154//              f2Act->setFacetNormal(f2Copy->getFacetNormal());
155                f2Act->setFacetNormal(f2Normal);
156                delete f2Normal;
157                f2Act->setUCN(f2Copy->getUCN());
158                f2Copy = f2Copy->next;
159        }       
160}
161               
162/** The default destructor */
163facet::~facet()
164{
165        if(this->fNormal!=NULL)
166                delete this->fNormal;
167        if(this->interiorPoint!=NULL)
168                delete this->interiorPoint;
169        /* Cleanup the codim2-structure */
170        if(this->codim==2)
171        {
172                facet *codim2Ptr;
173                codim2Ptr = this->codim2Ptr;
174                while(codim2Ptr!=NULL)
175                {
176                        if(codim2Ptr->fNormal!=NULL)
177                        {
178                                delete codim2Ptr->fNormal;
179                                codim2Ptr = codim2Ptr->next;
180                        }
181                }
182        }
183        if(this->codim2Ptr!=NULL)
184                delete this->codim2Ptr;
185        if(this->flipGB!=NULL)
186                idDelete((ideal *)&this->flipGB);
187        if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
188//              rDelete(this->flipRing);
189//      this->flipRing=NULL;
190        this->prev=NULL;
191        this->next=NULL;
192        //this=NULL;
193}
194               
195               
196/** \brief Comparison of facets*/
197inline bool gcone::areEqual(facet *f, facet *s)
198{
199        bool res = TRUE;
200        int notParallelCtr=0;
201        int ctr=0;
202        intvec* fNormal; //No new since ivCopy and therefore getFacetNormal return a new
203        intvec* sNormal;
204        fNormal = f->getFacetNormal();
205        sNormal = s->getFacetNormal();
206        if(!isParallel(fNormal,sNormal))
207                notParallelCtr++;
208        else//parallelity, so we check the codim2-facets
209        {
210                facet* f2Act;
211                facet* s2Act;
212                f2Act = f->codim2Ptr;           
213                ctr=0;
214                while(f2Act!=NULL)
215                {
216                        intvec* f2Normal;
217                        f2Normal = f2Act->getFacetNormal();
218                        s2Act = s->codim2Ptr;
219                        while(s2Act!=NULL)
220                        {
221                                intvec* s2Normal;
222                                s2Normal = s2Act->getFacetNormal();
223                                bool foo=areEqual(f2Normal,s2Normal);
224                                if(foo)
225                                        ctr++;
226                                s2Act = s2Act->next;
227                                delete s2Normal;
228                        }
229                        delete f2Normal;
230                        f2Act = f2Act->next;
231                }               
232        }
233        delete fNormal;
234        delete sNormal;
235        if(ctr==f->numCodim2Facets)
236                res=TRUE;
237        else
238                res=FALSE;                             
239        return res;
240}       
241               
242/** Stores the facet normal \param intvec*/
243inline void facet::setFacetNormal(intvec *iv)
244{
245        this->fNormal = ivCopy(iv);                     
246}
247               
248/** Hopefully returns the facet normal */
249inline intvec *facet::getFacetNormal()
250{                               
251        return ivCopy(this->fNormal);
252//      return this->fNormal;
253}
254               
255/** Method to print the facet normal*/
256inline void facet::printNormal()
257{
258        fNormal->show();
259}
260               
261/** Store the flipped GB*/
262inline void facet::setFlipGB(ideal I)
263{
264        this->flipGB=idCopy(I);
265}
266               
267/** Return the flipped GB
268Seems not be used
269Anyhow idCopy would make sense here.
270*/
271inline ideal facet::getFlipGB()
272{
273        return this->flipGB;
274}
275               
276/** Print the flipped GB*/
277inline void facet::printFlipGB()
278{
279#ifndef NDEBUG
280        idShow(this->flipGB);
281#endif
282}
283               
284/** Set the UCN */
285inline void facet::setUCN(int n)
286{
287        this->UCN=n;
288}
289               
290/** \brief Get the UCN
291* Returns the UCN iff this != NULL, else -1
292*/
293inline int facet::getUCN()
294{
295        if(this!=NULL)//  && this!=(facet *)0xfbfbfbfbfbfbfbfb)// || this!=(facet *)0xfbfbfbfb) )
296//      if(this!=NULL && ( this->fNormal!=(intvec *)0xfbfbfbfb || this->fNormal!=(intvec *)0xfbfbfbfbfbfbfbfb) )
297                return this->UCN;
298        else
299                return -1;
300}
301               
302/** Store an interior point of the facet */
303inline void facet::setInteriorPoint(intvec *iv)
304{
305        this->interiorPoint = ivCopy(iv);
306}
307               
308inline intvec *facet::getInteriorPoint()
309{
310        return this->interiorPoint;
311}       
312               
313/** \brief Debugging function
314* prints the facet normal an all (codim-2)-facets that belong to it
315*/
316inline void facet::fDebugPrint()
317{                       
318        facet *codim2Act;                       
319        codim2Act = this->codim2Ptr;
320        intvec *fNormal;
321        fNormal = this->getFacetNormal();       
322        cout << "=======================" << endl;
323        cout << "Facet normal = (";
324        fNormal->show(1,1);
325        cout << ")"<<endl;     
326        cout << "-----------------------" << endl;
327        cout << "Codim2 facets:" << endl;
328        while(codim2Act!=NULL)
329        {
330                intvec *f2Normal;
331                f2Normal = codim2Act->getFacetNormal();
332                cout << "(";
333                f2Normal->show(1,0);
334                cout << ")" << endl;
335                codim2Act = codim2Act->next;
336                delete f2Normal;
337        }
338        cout << "=======================" << endl;
339        delete fNormal; 
340}               
341               
342//friend class gcone;   //Bad style
343
344
345/**
346*\brief Implements the cone structure
347*
348* A cone is represented by a linked list of facet normals
349* @see facet
350*/
351
352
353/** \brief Default constructor.
354 *
355 * Initialises this->next=NULL and this->facetPtr=NULL
356 */
357gcone::gcone()
358{
359        this->next=NULL;
360        this->prev=NULL;
361        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
362        this->baseRing=currRing;
363        this->counter++;
364        this->UCN=this->counter;                       
365        this->numFacets=0;
366        this->ivIntPt=NULL;
367}
368               
369/** \brief Constructor with ring and ideal
370 *
371 * This constructor takes the root ring and the root ideal as parameters and stores
372 * them in the private members gcone::rootRing and gcone::inputIdeal
373 * This constructor is only called once in the computation of the Gröbner fan,
374 * namely for the very first cone. Therefore pred is set to 1.
375 * Might set it to this->UCN though...
376 * Since knowledge of the root ring is only needed when using reverse search,
377 * this constructor is not needed when using the "second" method
378*/
379gcone::gcone(ring r, ideal I)
380{
381        this->next=NULL;
382        this->prev=NULL;
383        this->facetPtr=NULL;                   
384//      this->rootRing=r;
385        this->inputIdeal=I;
386        this->baseRing=currRing;
387        this->counter++;
388        this->UCN=this->counter;
389        this->pred=1;
390        this->numFacets=0;
391        this->ivIntPt=NULL;
392}
393               
394/** \brief Copy constructor
395 *
396 * Copies a cone, sets this->gcBasis to the flipped GB
397 * Call this only after a successful call to gcone::flip which sets facet::flipGB
398*/             
399gcone::gcone(const gcone& gc, const facet &f)
400{
401        this->next=NULL;
402//      this->prev=(gcone *)&gc; //comment in to get a tree
403        this->prev=NULL;
404        this->numVars=gc.numVars;                                               
405        this->counter++;
406        this->UCN=this->counter;
407        this->pred=gc.UCN;
408        this->facetPtr=NULL;
409//      this->gcBasis=idCopy(f.flipGB);
410//      ring saveRing=currRing;
411//      ring tmpRing=f.flipRing;
412//      rComplete(f.flipRing);
413//      rChangeCurrRing(f.flipRing);
414        this->gcBasis=idrCopyR(f.flipGB, f.flipRing);
415//      this->inputIdeal=idCopy(this->gcBasis);
416        this->baseRing=rCopy(f.flipRing);
417//      rComplete(this->baseRing);
418//      rChangeCurrRing(saveRing);
419        this->numFacets=0;
420        this->ivIntPt=NULL;
421//      this->rootRing=NULL;
422        //rComplete(this->baseRing);
423        //rChangeCurrRing(this->baseRing);
424}
425               
426/** \brief Default destructor
427*/
428gcone::~gcone()
429{
430        if(this->gcBasis!=NULL)
431                idDelete((ideal *)&this->gcBasis);
432//      if(this->inputIdeal!=NULL)
433//              idDelete((ideal *)&this->inputIdeal);
434//      if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
435//              rDelete(this->rootRing);
436//      if(this->UCN!=1)
437//              rDelete(this->baseRing);
438        facet *fAct;
439        facet *fDel;
440        /*Delete the facet structure*/
441        fAct=this->facetPtr;
442        fDel=fAct;
443        while(fAct!=NULL)
444        {
445                fDel=fAct;
446                fAct=fAct->next;
447                delete fDel;
448        }
449        this->counter--;
450        //should be deleted in noRevS
451//      dd_FreeMatrix(this->ddFacets);
452        //dd_FreeMatrix(this->ddFacets);
453}                       
454       
455inline int gcone::getCounter()
456{
457        return this->counter;
458}
459       
460/** \brief Set the interior point of a cone */
461inline void gcone::setIntPoint(intvec *iv)
462{
463        this->ivIntPt=ivCopy(iv);
464}
465               
466/** \brief Return the interior point */
467inline intvec *gcone::getIntPoint()
468{
469        return ivCopy(this->ivIntPt);
470}
471               
472/** \brief Print the interior point */
473inline void gcone::showIntPoint()
474{
475        ivIntPt->show();
476}
477               
478/** \brief Print facets
479 * This is mainly for debugging purposes. Usually called from within gdb
480 */
481inline void gcone::showFacets(const short codim)
482{
483        facet *f;
484        switch(codim)
485        {
486                case 1:
487                        f = this->facetPtr;
488                        break;
489                case 2:
490                        f = this->facetPtr->codim2Ptr;
491                        break;
492        }       
493        while(f!=NULL)
494        {
495                intvec *iv;
496                iv = f->getFacetNormal();
497                cout << "(";
498                iv->show(1,0);                         
499                if(f->isFlippable==FALSE)
500                        cout << ")* ";
501                else
502                        cout << ") ";
503                delete iv;
504                f=f->next;
505        }
506        cout << endl;
507}
508               
509/** For debugging purposes only */
510inline volatile void gcone::showSLA(facet &f)
511{
512        facet *fAct;
513        fAct = &f;
514        facet *codim2Act;
515        codim2Act = fAct->codim2Ptr;
516       
517        cout << endl;
518        while(fAct!=NULL)
519        {
520                intvec *fNormal;               
521                fNormal=fAct->getFacetNormal();
522                cout << "(";
523                fNormal->show(1,0);
524                if(fAct->isFlippable==TRUE)
525                        cout << ") ";
526                else
527                        cout << ")* ";
528                delete fNormal;
529                codim2Act = fAct->codim2Ptr;
530                cout << " Codim2: ";
531                while(codim2Act!=NULL)
532                {
533                        intvec *f2Normal;
534                        f2Normal = codim2Act->getFacetNormal();
535                        cout << "(";
536                        f2Normal->show(1,0);
537                        cout << ") ";
538                        delete f2Normal;
539                        codim2Act = codim2Act->next;
540                }
541                cout << "UCN = " << fAct->getUCN() << endl;                             
542                fAct = fAct->next;
543        }
544}
545               
546inline void gcone::idDebugPrint(const ideal &I)
547{
548        int numElts=IDELEMS(I);
549        cout << "Ideal with " << numElts << " generators" << endl;
550        cout << "Leading terms: ";
551        for (int ii=0;ii<numElts;ii++)
552        {
553                pWrite0(pHead(I->m[ii]));
554                cout << ",";
555        }
556        cout << endl;
557}
558
559inline void gcone::invPrint(const ideal &I)
560{
561//      int numElts=IDELEMS(I);
562//      cout << "inv = ";
563//      for(int ii=0;ii<numElts;ii++);
564//      {
565//              pWrite0(pHead(I->m[ii]));
566//              cout << ",";
567//      }
568//      cout << endl;
569}
570
571inline bool gcone::isMonomial(const ideal &I)
572{
573        bool res = TRUE;
574        for(int ii=0;ii<IDELEMS(I);ii++)
575        {
576                if(pLength((poly)I->m[ii])>1)
577                {
578                        res = FALSE;
579                        break;
580                }                                               
581        }
582        return res;
583}
584               
585/** \brief Set gcone::numFacets */
586inline void gcone::setNumFacets()
587{
588}
589               
590/** \brief Get gcone::numFacets */
591inline int gcone::getNumFacets()
592{
593        return this->numFacets;
594}
595               
596inline int gcone::getUCN()
597{
598        if( this!=NULL)// && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
599                return this->UCN;
600        else
601                return -1;
602}
603
604inline int gcone::getPredUCN()
605{
606        return this->pred;
607}
608
609/** \brief Compute the normals of the cone
610 *
611 * This method computes a representation of the cone in terms of facet normals. It takes an ideal
612 * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
613 * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
614 * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
615 * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
616 * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
617 * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
618 *
619 * Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
620 * an interior point of the cone.
621                 */
622inline void gcone::getConeNormals(const ideal &I, bool compIntPoint)
623{
624#ifdef gfanp
625        timeval start, end;
626        gettimeofday(&start, 0);
627#endif
628        poly aktpoly;
629        int rows;                       // will contain the dimensions of the ineq matrix - deprecated by
630        dd_rowrange ddrows;
631        dd_colrange ddcols;
632        dd_rowset ddredrows;            // # of redundant rows in ddineq
633        dd_rowset ddlinset;             // the opposite
634        dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
635        dd_NumberType ddnumb=dd_Integer; //Number type
636        dd_ErrorType dderr=dd_NoError;         
637        //Compute the # inequalities i.e. rows of the matrix
638        rows=0; //Initialization
639        for (int ii=0;ii<IDELEMS(I);ii++)
640        {
641                aktpoly=(poly)I->m[ii];
642                rows=rows+pLength(aktpoly)-1;
643        }
644
645        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
646        ddrows=rows;
647        ddcols=this->numVars;
648        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
649        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
650               
651        // We loop through each g\in GB and compute the resulting inequalities
652        for (int i=0; i<IDELEMS(I); i++)
653        {
654                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
655                //simpler version of storing expvect diffs
656                int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
657//              int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
658                pGetExpV(aktpoly,leadexpv);
659                while(pNext(aktpoly)!=NULL)
660                {
661                        aktpoly=pNext(aktpoly);
662                        int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
663                        pGetExpV(aktpoly,tailexpv);
664                        for(int kk=1;kk<=this->numVars;kk++)
665                        {
666                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
667                        }
668                        aktmatrixrow += 1;
669                        omFree(tailexpv);
670                }               
671                omFree(leadexpv);       
672        } //for
673#if true
674        /*Let's make the preprocessing here. This could already be done in the above for-loop,
675        * but for a start it is more convenient here.
676        * We check the necessary condition of FJT p.18
677        * Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
678        */
679//      ideal initialForm=idInit(IDELEMS(I),1);
680        intvec *gamma=new intvec(this->numVars);
681        int falseGammaCounter=0;
682        int *redRowsArray=NULL;
683        int num_alloc=0;
684        int num_elts=0; 
685        for(int ii=0;ii<ddineq->rowsize;ii++)
686        {
687                ideal initialForm=idInit(IDELEMS(I),1);
688                //read row ii into gamma
689                double tmp;     
690                for(int jj=1;jj<=this->numVars;jj++)
691                {
692                        tmp=mpq_get_d(ddineq->matrix[ii][jj]);
693                        (*gamma)[jj-1]=(int)tmp;
694                }
695                computeInv((ideal&)I,initialForm,*gamma);
696                //Create leading ideal
697                ideal L=idInit(IDELEMS(initialForm),1);
698                for(int jj=0;jj<IDELEMS(initialForm);jj++)
699                {
700                        L->m[jj]=pCopy(pHead(initialForm->m[jj]));
701                }               
702               
703                LObject *P = new sLObject();
704                memset(P,0,sizeof(LObject));
705
706                for(int jj=0;jj<=IDELEMS(initialForm)-2;jj++)
707                {
708                        bool isMaybeFacet=FALSE;
709                        P->p1=initialForm->m[jj];       //build spolys of initialForm in_v
710
711                        for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
712                        {
713                                P->p2=initialForm->m[kk];
714                                ksCreateSpoly(P);                                                       
715                                if(P->p!=NULL)  //spoly non zero=?
716                                {       
717                                        poly p=pInit();                 
718                                        poly q=pInit();
719                                        p=pCopy(P->p);
720                                        q=pHead(p);     //Monomial q
721                                        isMaybeFacet=FALSE;
722                                        while(p!=NULL)
723                                        {
724                                                q=pHead(p);                                             
725//                                              unsigned long sevSpoly=pGetShortExpVector(q);
726//                                              unsigned long not_sevL;                                         
727                                                for(int ll=0;ll<IDELEMS(L);ll++)
728                                                {
729//                                                      not_sevL=~pGetShortExpVector(L->m[ll]);//                                       
730                                                        //if(!(sevSpoly & not_sevL) && pLmDivisibleBy(L->m[ll],q) )//i.e. spoly is in L
731                                                        if(pLmEqual(L->m[ll],q) || pDivisibleBy(L->m[ll],q))
732                                                        {                                                       
733                                                                isMaybeFacet=TRUE;
734                                                                break;//for
735                                                        }
736                                                }
737                                                if(isMaybeFacet==TRUE)
738                                                {                                                       
739                                                        break;//while(p!=NULL)
740                                                }
741                                                p=pNext(p);                                             
742                                        }//while
743                                        pDelete(&p);
744                                        pDelete(&q);
745                                        if(isMaybeFacet==FALSE)
746                                        {
747                                                dd_set_si(ddineq->matrix[ii][0],1);                                             
748                                                if(num_alloc==0)
749                                                        num_alloc += 1;
750                                                else
751                                                        num_alloc += 1;
752                                                void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
753                                                if(!tmp)
754                                                {
755                                                        WerrorS("Woah dude! Couldn't realloc memory\n");
756                                                        exit(-1);
757                                                }
758                                                redRowsArray = (int*)tmp;
759                                                redRowsArray[num_elts]=ii;
760                                                num_elts++;
761                                                //break;//for(int kk, since we have found one that is not in L 
762                                                goto _start;    //mea culpa, mea culpa, mea maxima culpa
763                                        }
764                                }//if(P->p!=NULL)                               
765                        }//for k
766                }//for jj
767                _start:;
768                idDelete(&L);
769                delete P;
770                idDelete(&initialForm);
771                //idDelete(L);         
772        }//for(ii<ddineq-rowsize
773        delete gamma;
774        int offset=0;//needed for correction of redRowsArray[ii]
775        for( int ii=0;ii<num_elts;ii++ )
776        {
777                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
778                offset++;
779        }
780        free(redRowsArray);
781        /*And now for the strictly positive rows
782        * Doesn't gain significant speedup
783        */
784        /*int *posRowsArray=NULL;
785        num_alloc=0;
786        num_elts=0;
787        for(int ii=0;ii<ddineq->rowsize;ii++)
788        {
789                intvec *ivPos = new intvec(this->numVars);
790                for(int jj=0;jj<this->numVars;jj++)
791                        (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
792                bool isStrictlyPos=FALSE;
793                int posCtr=0;           
794                for(int jj=0;jj<this->numVars;jj++)
795                {
796                        intvec *ivCanonical = new intvec(this->numVars);
797                        jj==0 ? (*ivCanonical)[ivPos->length()-1]=1 : (*ivCanonical)[jj-1]=1;
798                        if(dotProduct(*ivCanonical,*ivPos)!=0)
799                        {
800                                if ((*ivPos)[jj]>=0)
801                                {                               
802                                        posCtr++;                               
803                                }
804                        }                       
805                        delete ivCanonical;
806                }
807                if(posCtr==ivPos->length())
808                        isStrictlyPos=TRUE;
809                if(isStrictlyPos==TRUE)
810                {
811                        if(num_alloc==0)
812                                num_alloc += 1;
813                        else
814                                num_alloc += 1;
815                        void *tmp = realloc(posRowsArray,(num_alloc*sizeof(int)));
816                        if(!tmp)
817                        {
818                                WerrorS("Woah dude! Couldn't realloc memory\n");
819                                exit(-1);
820                        }
821                        posRowsArray = (int*)tmp;
822                        posRowsArray[num_elts]=ii;
823                        num_elts++;     
824                }
825                delete ivPos;
826        }
827        offset=0;
828        for(int ii=0;ii<num_elts;ii++)
829        {
830                dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
831                offset++;
832        }
833        free(posRowsArray);*/
834#endif
835
836        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
837        //time(&canonicalizeTac);
838        //cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl;
839        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
840        ddcols = ddineq->colsize;
841                       
842        //ddCreateMatrix(ddrows,ddcols+1);
843        ddFacets = dd_CopyMatrix(ddineq);
844                       
845        /*Write the normals into class facet*/ 
846        facet *fAct;    //pointer to active facet       
847        int numNonFlip=0;
848        for (int kk = 0; kk<ddrows; kk++)
849        {
850                intvec *load = new intvec(this->numVars);       //intvec to store a single facet normal that will then be stored via setFacetNormal
851                for (int jj = 1; jj <ddcols; jj++)
852                {
853                        double foo;
854                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
855                        (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load           
856                }//for (int jj = 1; jj <ddcols; jj++)
857                               
858                /*Quick'n'dirty hack for flippability*/ 
859                bool isFlip=FALSE;                                                             
860                for (int jj = 0; jj<load->length(); jj++)
861                {
862                        intvec *ivCanonical = new intvec(load->length());
863                        (*ivCanonical)[jj]=1;
864//                                      cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl;
865                        if (dotProduct(*load,*ivCanonical)<0)   
866                                        //if (ivMult(load,ivCanonical)<0)
867                        {
868                                isFlip=TRUE;
869                                break;  //URGHS
870                        }
871                        delete ivCanonical;
872                }/*End of check for flippability*/
873                if (isFlip==FALSE)
874                {
875                        this->numFacets++;
876                        numNonFlip++;
877                        if(this->numFacets==1)
878                        {
879                                facet *fRoot = new facet();
880                                this->facetPtr = fRoot;
881                                fAct = fRoot;
882                                               
883                        }
884                        else
885                        {
886                                fAct->next = new facet();
887                                fAct = fAct->next;
888                        }
889                        fAct->isFlippable=FALSE;
890                        fAct->setFacetNormal(load);
891                        fAct->setUCN(this->getUCN());
892#ifdef gfan_DEBUG
893                        cout << "Marking facet (";
894                        load->show(1,0);
895                        cout << ") as non flippable" << endl;                           
896#endif
897                }
898                else
899                {       /*Now load should be full and we can call setFacetNormal*/                                     
900                        this->numFacets++;
901                        if(this->numFacets==1)
902                        {
903                                facet *fRoot = new facet();
904                                this->facetPtr = fRoot;
905                                fAct = fRoot;
906                        }
907                        else
908                        {
909                                fAct->next = new facet();
910                                fAct = fAct->next;
911                        }
912                        fAct->isFlippable=TRUE;
913                        fAct->setFacetNormal(load);
914                        fAct->setUCN(this->getUCN());                                   
915                }//if (isFlippable==FALSE)
916                delete load;                           
917        }//for (int kk = 0; kk<ddrows; kk++)
918                       
919        //In cases like I=<x-1,y-1> there are only non-flippable facets...
920        if(numNonFlip==this->numFacets)
921        {                                       
922                WerrorS ("Only non-flippable facets. Terminating...\n");
923//              exit(-1);//Bit harsh maybe...
924        }
925                       
926        /*
927        Now we should have a linked list containing the facet normals of those facets that are
928        -irredundant
929        -flipable
930        Adressing is done via *facetPtr
931        */                     
932        if (compIntPoint==TRUE)
933        {
934                intvec *iv = new intvec(this->numVars);
935                dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
936                int jj=1;
937                for (int ii=0;ii<=this->numVars;ii++)
938                {
939                        dd_set_si(posRestr->matrix[ii][jj],1);
940                        jj++;                                                   
941                }
942                dd_MatrixAppendTo(&ddineq,posRestr);
943                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
944                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
945                delete iv;
946                dd_FreeMatrix(posRestr);
947        }       
948        //Clean up but don't delete the return value!   
949        dd_FreeMatrix(ddineq);
950        set_free(ddredrows);
951        set_free(ddlinset);
952        free(ddnewpos);
953#ifdef gfanp
954        gettimeofday(&end, 0);
955        time_getConeNormals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
956#endif
957       
958
959}//gcone::getConeNormals(ideal I)
960               
961/** \brief Compute the (codim-2)-facets of a given cone
962 * This method is used during noRevS
963 * Additionally we check whether the codim2-facet normal is strictly positive. Otherwise
964 * the facet is marked as non-flippable.
965 */
966inline void gcone::getCodim2Normals(const gcone &gc)
967{
968#ifdef gfanp
969        timeval start, end;
970        gettimeofday(&start, 0);
971#endif
972        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
973        facet *fAct;
974        fAct = this->facetPtr;         
975        facet *codim2Act;
976        //codim2Act = this->facetPtr->codim2Ptr;
977        dd_MatrixPtr ddineq;//,P,ddakt;
978        dd_ErrorType err;
979        //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
980        ddineq = dd_CopyMatrix(gc.ddFacets);
981        /*Now set appropriate linearity*/
982        for (int ii=0; ii<this->numFacets; ii++)                       
983        {       
984//              dd_rowset LL;
985//              set_initialize(&LL,ddineq->rowsize);
986//              set_copy(LL,ddineq->linset);
987//              set_copy(ddineq->linset,LL);
988                dd_rowset impl_linset, redset;
989                dd_rowindex newpos;
990                dd_MatrixPtr ddakt;
991                ddakt = dd_CopyMatrix(ddineq);
992                ddakt->representation=dd_Inequality;
993                set_addelem(ddakt->linset,ii+1);/*Now set appropriate linearity*/
994#ifdef gfanp
995                timeval t_ddMC_start, t_ddMC_end;
996                gettimeofday(&t_ddMC_start,0);
997#endif                         
998                /*dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);*/
999//              set_copy(LL,ddakt->linset);
1000                dd_PolyhedraPtr ddpolyh;
1001                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
1002                dd_MatrixPtr P;
1003                P=dd_CopyGenerators(ddpolyh);
1004                dd_FreePolyhedra(ddpolyh);
1005#ifdef gfanp
1006                gettimeofday(&t_ddMC_end,0);
1007                t_ddMC += (t_ddMC_end.tv_sec - t_ddMC_start.tv_sec + 1e-6*(t_ddMC_end.tv_usec - t_ddMC_start.tv_usec));
1008#endif 
1009                /* We loop through each row of P normalize it by making all
1010                * entries integer ones and add the resulting vector to the
1011                * int matrix facet::codim2Facets */
1012                for (int jj=1;jj<=/*ddakt*/P->rowsize;jj++)
1013                {                                       
1014                        fAct->numCodim2Facets++;
1015                        if(fAct->numCodim2Facets==1)                                   
1016                        {                                               
1017                                fAct->codim2Ptr = new facet(2);                                         
1018                                codim2Act = fAct->codim2Ptr;
1019                        }
1020                        else
1021                        {
1022                                codim2Act->next = new facet(2);
1023                                codim2Act = codim2Act->next;
1024                        }
1025                        intvec *n = new intvec(this->numVars);
1026#ifdef gfanp
1027                        timeval t_mI_start, t_mI_end;
1028                        gettimeofday(&t_mI_start,0);
1029#endif
1030                        makeInt(P,jj,*n);
1031                        /*for(int kk=0;kk<this->numVars;kk++)
1032                        {
1033                                int foo;
1034                                foo = (int)mpq_get_d(ddakt->matrix[ii][kk+1]);
1035                                (*n)[kk]=foo;
1036                        }*/
1037#ifdef gfanp
1038                        gettimeofday(&t_mI_end,0);
1039                        t_mI += (t_mI_end.tv_sec - t_mI_start.tv_sec + 1e-6*(t_mI_end.tv_usec - t_mI_start.tv_usec));
1040#endif
1041                        codim2Act->setFacetNormal(n);
1042                        delete n;                                       
1043                }               
1044                /*We check whether the facet spanned by the codim-2 facets
1045                * intersects with the positive orthant. Otherwise we define this
1046                * facet to be non-flippable. Works since we set the appropriate
1047                * linearity for ddakt above.
1048                */
1049                intvec *iv_intPoint = new intvec(this->numVars);
1050                dd_MatrixPtr shiftMatrix;
1051                dd_MatrixPtr intPointMatrix;
1052                shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
1053                for(int kk=0;kk<this->numVars;kk++)
1054                {
1055                        dd_set_si(shiftMatrix->matrix[kk][0],1);
1056                        dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
1057                }
1058                intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
1059#ifdef gfanp
1060                timeval t_iP_start, t_iP_end;
1061                gettimeofday(&t_iP_start, 0);
1062#endif         
1063                interiorPoint(intPointMatrix,*iv_intPoint);
1064//              dd_rowset impl_linste,lbasis;
1065//              dd_LPSolutionPtr lps=NULL;
1066//              dd_ErrorType err;
1067//              dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
1068#ifdef gfanp
1069                gettimeofday(&t_iP_end, 0);
1070                t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
1071#endif
1072                for(int ll=0;ll<this->numVars;ll++)
1073                {
1074                        if( (*iv_intPoint)[ll] < 0 )
1075                        {
1076                                fAct->isFlippable=FALSE;
1077                                break;
1078                        }
1079                }
1080                /*End of check*/                                       
1081                fAct = fAct->next;     
1082                dd_FreeMatrix(ddakt);
1083                dd_FreeMatrix(shiftMatrix);
1084                dd_FreeMatrix(intPointMatrix);         
1085                delete iv_intPoint;
1086                dd_FreeMatrix(P);
1087//              set_free(impl_linset);
1088//              set_free(redset);               
1089//              free(newpos);
1090//              set_free(LL);
1091        }//for 
1092        dd_FreeMatrix(ddineq);
1093#ifdef gfanp
1094        gettimeofday(&end, 0);
1095        time_getCodim2Normals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1096#endif
1097}
1098               
1099/** \brief Compute the Groebner Basis on the other side of a shared facet
1100 *
1101 * Implements algorithm 4.3.2 from Anders' thesis.
1102 * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
1103 * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
1104 * is parallel to \f$ leadexp(g) \f$
1105 * Parallelity is checked using basic linear algebra. See gcone::isParallel.
1106 * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
1107 * computing an interior point of the facet and taking all terms having the same weight with respect
1108 * to this interior point.
1109 *\param ideal, facet
1110 * Input: a marked,reduced Groebner basis and a facet
1111 */
1112inline void gcone::flip(ideal gb, facet *f)             //Compute "the other side"
1113{       
1114#ifdef gfanp
1115        timeval start, end;
1116        gettimeofday(&start, 0);
1117#endif         
1118        intvec *fNormal;// = new intvec(this->numVars); //facet normal, check for parallelity                   
1119        fNormal = f->getFacetNormal();  //read this->fNormal;
1120
1121//      std::cout << "running gcone::flip" << std::endl;
1122//      std::cout << "flipping UCN " << this->getUCN() << endl;
1123//      cout << "over facet (";
1124//      fNormal->show(1,0);
1125//      cout << ") with UCN " << f->getUCN();
1126//      std::cout << std::endl;
1127        if(this->getUCN() != f->getUCN())
1128        {
1129                WerrorS("Uh oh... Trying to flip over facet with incompatible UCN");
1130                exit(-1);
1131        }
1132        /*1st step: Compute the initial ideal*/
1133        /*poly initialFormElement[IDELEMS(gb)];*/       //array of #polys in GB to store initial form
1134        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
1135       
1136        computeInv(gb,initialForm,*fNormal);
1137
1138#ifdef gfan_DEBUG
1139/*      cout << "Initial ideal is: " << endl;
1140        idShow(initialForm);
1141        //f->printFlipGB();*/
1142//      cout << "===" << endl;
1143#endif                 
1144        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
1145        /*Substep 2.1
1146        compute $G_{-\alpha}(in_v(I))
1147        see journal p. 66
1148        NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
1149        srcRing already has a weighted ordering
1150        */
1151        ring srcRing=currRing;
1152        ring tmpRing;
1153                       
1154        if( (srcRing->order[0]!=ringorder_a))
1155        {
1156                intvec *iv = new intvec(this->numVars);
1157                iv = ivNeg(fNormal);
1158//              tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
1159                tmpRing=rCopyAndAddWeight(srcRing,iv);
1160                delete iv;
1161        }
1162        else
1163        {
1164                tmpRing=rCopy0(srcRing);
1165                int length=fNormal->length();
1166                int *A=(int *)omAlloc0(length*sizeof(int));
1167                for(int jj=0;jj<length;jj++)
1168                {
1169                        A[jj]=-(*fNormal)[jj];
1170                }
1171                omFree(tmpRing->wvhdl[0]);
1172                tmpRing->wvhdl[0]=(int*)A;
1173                tmpRing->block1[0]=length;
1174                rComplete(tmpRing);
1175                //omFree(A);
1176        }
1177        delete fNormal; 
1178        rChangeCurrRing(tmpRing);       
1179                       
1180        ideal ina;                     
1181        ina=idrCopyR(initialForm,srcRing);
1182        idDelete(&initialForm);
1183        ideal H;
1184//      H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
1185        H=kStd(ina,NULL,testHomog,NULL);        //This is \mathcal(G)_{>_-\alpha}(in_v(I))
1186        idSkipZeroes(H);
1187        idDelete(&ina);
1188
1189        /*Substep 2.2
1190        do the lifting and mark according to H
1191        */
1192        rChangeCurrRing(srcRing);
1193        ideal srcRing_H;
1194        ideal srcRing_HH;                       
1195        srcRing_H=idrCopyR(H,tmpRing);
1196        //H is needed further below, so don't idDelete here
1197        srcRing_HH=ffG(srcRing_H,this->gcBasis);
1198        idDelete(&srcRing_H);
1199       
1200        /*Substep 2.2.1
1201         * Mark according to G_-\alpha
1202         * Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
1203         * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
1204         * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
1205         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
1206         * compute the difference accordingly
1207        */
1208#ifdef gfanp
1209        timeval t_markings_start, t_markings_end;
1210        gettimeofday(&t_markings_start, 0);
1211#endif         
1212        bool markingsAreCorrect=FALSE;
1213        dd_MatrixPtr intPointMatrix;
1214        int iPMatrixRows=0;
1215        dd_rowrange aktrow=0;                   
1216        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1217        {
1218                poly aktpoly=(poly)srcRing_HH->m[ii];//This is a pointer, so don't pDelete
1219                iPMatrixRows = iPMatrixRows+pLength(aktpoly);           
1220        }
1221        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
1222         * construction of the differences
1223         */
1224        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); 
1225        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
1226       
1227        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1228        {
1229                markingsAreCorrect=FALSE;       //crucial to initialise here
1230                poly aktpoly=srcRing_HH->m[ii]; //Only a pointer, so don't pDelete
1231                /*Comparison of leading monomials is done via exponent vectors*/
1232                for (int jj=0;jj<IDELEMS(H);jj++)
1233                {
1234                        int *src_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
1235                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
1236                        pGetExpV(aktpoly,src_ExpV);
1237                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
1238                        pGetExpV(pCopy(H->m[ii]),dst_ExpV);
1239                        rChangeCurrRing(srcRing);
1240                        bool expVAreEqual=TRUE;
1241                        for (int kk=1;kk<=this->numVars;kk++)
1242                        {
1243#ifdef gfan_DEBUG
1244//                                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
1245#endif
1246                                if (src_ExpV[kk]!=dst_ExpV[kk])
1247                                {
1248                                        expVAreEqual=FALSE;
1249                                }
1250                        }                                       
1251                                        //if (*src_ExpV == *dst_ExpV)
1252                        if (expVAreEqual==TRUE)
1253                        {
1254                                markingsAreCorrect=TRUE; //everything is fine
1255#ifdef gfan_DEBUG
1256//                              cout << "correct markings" << endl;
1257#endif
1258                        }//if (pHead(aktpoly)==pHead(H->m[jj])
1259                        omFree(src_ExpV);
1260                        omFree(dst_ExpV);
1261                }//for (int jj=0;jj<IDELEMS(H);jj++)
1262               
1263                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
1264                if (markingsAreCorrect==TRUE)
1265                {
1266                        pGetExpV(aktpoly,leadExpV);
1267                }
1268                else
1269                {
1270                        rChangeCurrRing(tmpRing);
1271                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
1272                        rChangeCurrRing(srcRing);
1273                }
1274                /*compute differences of the expvects*/                         
1275                while (pNext(aktpoly)!=NULL)
1276                {
1277                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
1278                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
1279                        is not omitted when computing the differences*/
1280                        if(markingsAreCorrect==TRUE)
1281                        {
1282                                aktpoly=pNext(aktpoly);
1283                                pGetExpV(aktpoly,v);
1284                        }
1285                        else
1286                        {
1287                                pGetExpV(pHead(aktpoly),v);
1288                                markingsAreCorrect=TRUE;
1289                        }
1290
1291                        for (int jj=0;jj<this->numVars;jj++)
1292                        {                               
1293                                /*Store into ddMatrix*/                                                         
1294                                dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
1295                        }
1296                        aktrow +=1;
1297                        omFree(v);
1298                }
1299//              omFree(v);
1300                omFree(leadExpV);
1301        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1302#ifdef gfanp
1303        gettimeofday(&t_markings_end, 0);
1304        t_markings += (t_markings_end.tv_sec - t_markings_start.tv_sec + 1e-6*(t_markings_end.tv_usec - t_markings_start.tv_usec));
1305#endif
1306        /*Now it is safe to idDelete(H)*/
1307        idDelete(&H);
1308        /*Now we add the constraint for the standard simplex*/
1309// #ifdef gfanp
1310//      timeval t_dd_start, t_dd_end;
1311//      gettimeofday(&t_dd_start, 0);
1312// #endif
1313        dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
1314        for (int jj=1;jj<=this->numVars;jj++)
1315        {
1316                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
1317        }
1318        //Let's make sure we compute interior points from the positive orthant
1319        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1320       
1321        int jj=1;
1322        for (int ii=0;ii<this->numVars;ii++)
1323        {
1324                dd_set_si(posRestr->matrix[ii][jj],1);
1325                jj++;                                                   
1326        }
1327        dd_MatrixAppendTo(&intPointMatrix,posRestr);
1328        dd_FreeMatrix(posRestr);
1329        /*Insert preprocessing here. If it is before the dd_CreateMatrix things go pear shaped
1330        * Otherwise things will go pear shaped if called after the standard simplex constraints are added
1331        */
1332        preprocessInequalities(intPointMatrix);
1333        intvec *iv_weight = new intvec(this->numVars);
1334#ifdef gfanp
1335        timeval t_dd_start, t_dd_end;
1336        gettimeofday(&t_dd_start, 0);
1337#endif
1338        dd_ErrorType err;
1339        dd_rowset implLin, redrows;
1340        dd_rowindex newpos;
1341       
1342        dd_MatrixCanonicalize(&intPointMatrix,&implLin,&redrows,&newpos,&err);
1343//      dd_MatrixCanonicalizeLinearity(&intPointMatrix,&implLin, &newpos,&err);
1344        //dd_MatrixRedundancyRemove is our time sink!
1345//      dd_MatrixRedundancyRemove(&intPointMatrix,&redrows,&newpos,&err);
1346        interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
1347        dd_FreeMatrix(intPointMatrix);
1348        /*Crude attempt for interior point */
1349        /*dd_PolyhedraPtr ddpolyh;
1350        dd_ErrorType err;
1351        dd_rowset impl_linset,redset;
1352        dd_rowindex newpos;
1353        dd_MatrixCanonicalize(&intPointMatrix, &impl_linset, &redset, &newpos, &err);
1354        ddpolyh=dd_DDMatrix2Poly(intPointMatrix, &err);
1355        dd_MatrixPtr P;
1356        P=dd_CopyGenerators(ddpolyh);
1357        dd_FreePolyhedra(ddpolyh);
1358        for(int ii=0;ii<P->rowsize;ii++)
1359        {
1360                intvec *iv_row=new intvec(this->numVars);
1361                makeInt(P,ii+1,*iv_row);
1362                iv_weight =ivAdd(iv_weight, iv_row);
1363                delete iv_row;
1364        }
1365        dd_FreeMatrix(P);
1366        dd_FreeMatrix(intPointMatrix);*/
1367#ifdef gfanp
1368        gettimeofday(&t_dd_end, 0);
1369        t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec));
1370#endif                 
1371        /*Step 3
1372        turn the minimal basis into a reduced one
1373        */                     
1374        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
1375        // Thus:
1376        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
1377        ring dstRing=rCopy0(tmpRing);
1378        int length=iv_weight->length();
1379        int *A=(int *)omAlloc0(length*sizeof(int));
1380        for(int jj=0;jj<length;jj++)
1381        {
1382                A[jj]=(*iv_weight)[jj];
1383        }
1384        dstRing->wvhdl[0]=(int*)A;
1385        rComplete(dstRing);
1386//      omFree(A);                                     
1387        rChangeCurrRing(dstRing);
1388        rDelete(tmpRing);
1389        delete iv_weight;
1390
1391        ideal dstRing_I;                       
1392        dstRing_I=idrCopyR(srcRing_HH,srcRing);
1393        idDelete(&srcRing_HH); //Hmm.... causes trouble - no more
1394        //dstRing_I=idrCopyR(inputIdeal,srcRing);
1395        BITSET save=test;
1396        test|=Sy_bit(OPT_REDSB);
1397        test|=Sy_bit(OPT_REDTAIL);
1398#ifdef gfan_DEBUG
1399        test|=Sy_bit(6);        //OPT_DEBUG
1400#endif
1401        ideal tmpI;
1402        //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
1403        //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
1404        tmpI = dstRing_I;
1405        dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
1406        idDelete(&tmpI);
1407        idNorm(dstRing_I);                     
1408//      kInterRed(dstRing_I);
1409        idSkipZeroes(dstRing_I);
1410        test=save;
1411        /*End of step 3 - reduction*/
1412                       
1413        f->setFlipGB(dstRing_I);//store the flipped GB
1414        idDelete(&dstRing_I);
1415        f->flipRing=rCopy(dstRing);     //store the ring on the other side
1416//#ifdef gfan_DEBUG
1417//      cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
1418//      this->idDebugPrint(dstRing_I);
1419//      cout << endl;
1420//#endif                       
1421        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
1422        rDelete(dstRing);
1423#ifdef gfanp
1424        gettimeofday(&end, 0);
1425        time_flip += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1426#endif
1427}//void flip(ideal gb, facet *f)
1428
1429/** \brief Compute initial ideal
1430 * Compute the initial ideal in_v(G) wrt a (possible) facet normal
1431 * used in gcone::getFacetNormal in order to preprocess possible facet normals
1432 * and in gcone::flip for obvious reasons.
1433*/
1434inline void gcone::computeInv(ideal &gb, ideal &initialForm, intvec &fNormal)
1435{
1436#ifdef gfanp
1437        timeval start, end;
1438        gettimeofday(&start, 0);
1439#endif
1440        for (int ii=0;ii<IDELEMS(gb);ii++)
1441        {
1442                poly initialFormElement;
1443                poly aktpoly = (poly)gb->m[ii];//Ptr, so don't pDelete(aktpoly)
1444                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
1445                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
1446//              initialFormElement[ii]=pHead(aktpoly);
1447                initialFormElement=pHead(aktpoly);
1448                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
1449                {
1450                        intvec *check = new intvec(this->numVars);
1451                        aktpoly=pNext(aktpoly); //next term
1452//                      pSetm(aktpoly);
1453                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
1454                        pGetExpV(aktpoly,v);           
1455                        /* Convert (int)v into (intvec)check */                 
1456                        for (int jj=0;jj<this->numVars;jj++)
1457                        {
1458                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
1459                        }
1460                        if (isParallel(*check,fNormal)) //pass *check when
1461                        {
1462                                //Found a parallel vector. Add it
1463//                              initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
1464                                initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));
1465                        }
1466                        omFree(v);
1467                        delete check;
1468                }//while
1469#ifdef gfan_DEBUG
1470//              cout << "Initial Form=";                               
1471//              pWrite(initialFormElement[ii]);
1472//              cout << "---" << endl;
1473#endif
1474                /*Now initialFormElement must be added to (ideal)initialForm */
1475//              initialForm->m[ii]=pCopy(initialFormElement[ii]);
1476//              pDelete(&initialFormElement[ii]);
1477                initialForm->m[ii]=pCopy(initialFormElement);
1478                pDelete(&initialFormElement);
1479                omFree(leadExpV);
1480//              delete check;
1481        }//for
1482#ifdef gfanp
1483        gettimeofday(&end, 0);
1484        time_computeInv += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1485#endif
1486}
1487
1488/** \brief Compute the remainder of a polynomial by a given ideal
1489 *
1490 * Compute \f$ f^{\mathcal{G}} \f$
1491 * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
1492 * However, since we are only interested in the remainder, there is no need to
1493 * compute the factors \f$ a_i \f$
1494 */
1495//NOTE: Should be replaced by kNF or kNF2
1496//NOTE: Done
1497//NOTE: removed with r12286
1498               
1499/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
1500*/
1501//NOTE: use kNF or kNF2 instead of restOfDivision
1502inline ideal gcone::ffG(const ideal &H, const ideal &G)
1503{
1504//                      cout << "Entering ffG" << endl;
1505        int size=IDELEMS(H);
1506        ideal res=idInit(size,1);
1507//      poly temp1;//=pInit();
1508//      poly temp2;//=pInit();
1509//      poly temp3;//=pInit();  //polys to temporarily store values for pSub
1510        for (int ii=0;ii<size;ii++)
1511        {
1512                poly temp1=pInit();
1513                poly temp2=pInit();
1514                poly temp3=pInit();
1515//              res->m[ii]=restOfDiv(H->m[ii],G);
1516//              res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
1517                temp1=pCopy(H->m[ii]);
1518//              temp2=pCopy(res->m[ii]);
1519                //NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
1520                temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));
1521                temp3=pSub(temp1, temp2);
1522                res->m[ii]=pCopy(temp3);
1523                //res->m[ii]=pSub(temp1,temp2); //buggy         
1524                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
1525                pDelete(&temp1);
1526//              pDelete(&temp2);
1527//              pDelete(&temp3); //NOTE does not work, so commented out
1528        }
1529        return res;
1530}
1531       
1532/** \brief Preprocessing of inequlities
1533* Do some preprocessing on the matrix of inequalities
1534* 1) Replace several constraints on the pos. orthants by just one for each orthant
1535* 2) Remove duplicates of inequalities
1536* 3) Remove inequalities that arise as sums of other inequalities
1537*/
1538void gcone::preprocessInequalities(dd_MatrixPtr &ddineq)
1539{
1540//Remove strictly positive rows
1541//      int *posRowsArray=NULL;
1542//      int num_alloc=0;
1543//      int num_elts=0;
1544//      for(int ii=0;ii<ddineq->rowsize;ii++)
1545//      {
1546//              intvec *ivPos = new intvec(pVariables);
1547//              for(int jj=0;jj<pVariables;jj++)
1548//                      (*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
1549//              bool isStrictlyPos=FALSE;
1550//              int posCtr=0;           
1551//              for(int jj=0;jj<pVariables;jj++)
1552//              {
1553//                      intvec *ivCanonical = new intvec(pVariables);
1554//                      jj==0 ? (*ivCanonical)[ivPos->length()-1]=1 : (*ivCanonical)[jj-1]=1;
1555//                      if(dotProduct(*ivCanonical,*ivPos)!=0)
1556//                      {
1557//                              if ((*ivPos)[jj]>=0)
1558//                              {                               
1559//                                      posCtr++;                               
1560//                              }
1561//                      }                       
1562//                      delete ivCanonical;
1563//              }
1564//              if(posCtr==ivPos->length())
1565//                      isStrictlyPos=TRUE;
1566//              if(isStrictlyPos==TRUE)
1567//              {
1568//                      if(num_alloc==0)
1569//                              num_alloc += 1;
1570//                      else
1571//                              num_alloc += 1;
1572//                      void *tmp = realloc(posRowsArray,(num_alloc*sizeof(int)));
1573//                      if(!tmp)
1574//                      {
1575//                              WerrorS("Woah dude! Couldn't realloc memory\n");
1576//                              exit(-1);
1577//                      }
1578//                      posRowsArray = (int*)tmp;
1579//                      posRowsArray[num_elts]=ii;
1580//                      num_elts++;     
1581//              }
1582//              delete ivPos;
1583//      }
1584        int offset=0;
1585//      for(int ii=0;ii<num_elts;ii++)
1586//      {
1587//              dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
1588//              offset++;
1589//      }
1590//      free(posRowsArray);
1591        //Remove zeroes
1592        int rowsize=ddineq->rowsize;
1593        for(int ii=0;ii<ddineq->rowsize;ii++)
1594        {
1595                intvec *iv = new intvec(this->numVars);
1596                int posCtr=0;
1597                for(int jj=0;jj<this->numVars;jj++)
1598                {
1599                        (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
1600                        if((*iv)[ii]>0)
1601                                posCtr++;
1602                }
1603                if( (iv->compare(0)==0)) //|| (posCtr==iv->length()) )
1604                {
1605                        dd_MatrixRowRemove(&ddineq,ii+1);
1606                        ii--;           
1607//                      rowsize=ddineq->rowsize;
1608                }
1609                delete iv;
1610        }
1611        //Remove duplicates of rows
1612//      posRowsArray=NULL;
1613//      num_alloc=0;
1614//      num_elts=0;
1615//      offset=0;
1616//      int num_newRows = ddineq->rowsize;
1617//      for(int ii=0;ii<ddineq->rowsize-1;ii++)
1618//      for(int ii=0;ii<num_newRows-1;ii++)
1619//      {
1620//              intvec *iv = new intvec(this->numVars);//1st vector to check against
1621//              for(int jj=0;jj<this->numVars;jj++)
1622//                      (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
1623//              for(int jj=ii+1;jj</*ddineq->rowsize*/num_newRows;jj++)
1624//              {
1625//                      intvec *ivCheck = new intvec(this->numVars);//Checked against iv
1626//                      for(int kk=0;kk<this->numVars;kk++)
1627//                              (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
1628//                      if (iv->compare(ivCheck)==0)
1629//                      {
1630// //                           cout << "=" << endl;
1631// //                           num_alloc++;
1632// //                           void *tmp=realloc(posRowsArray,(num_alloc*sizeof(int)));
1633// //                           if(!tmp)
1634// //                           {
1635// //                                   WerrorS("Woah dude! Couldn't realloc memory\n");
1636// //                                   exit(-1);
1637// //                           }
1638// //                           posRowsArray = (int*)tmp;
1639// //                           posRowsArray[num_elts]=jj;
1640// //                           num_elts++;
1641//                              dd_MatrixRowRemove(&ddineq,jj+1);
1642//                              num_newRows = ddineq->rowsize;
1643//                      }
1644//                      delete ivCheck;
1645//              }
1646//              delete iv;
1647//      }
1648//      for(int ii=0;ii<num_elts;ii++)
1649//      {
1650//              dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
1651//              offset++;
1652//      }
1653//      free(posRowsArray);
1654        //Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/       
1655}//preprocessInequalities
1656       
1657/** \brief Compute a Groebner Basis
1658 *
1659 * Computes the Groebner basis and stores the result in gcone::gcBasis
1660 *\param ideal
1661 *\return void
1662 */
1663inline void gcone::getGB(const ideal &inputIdeal)               
1664{
1665        BITSET save=test;
1666        test|=Sy_bit(OPT_REDSB);
1667        test|=Sy_bit(OPT_REDTAIL);
1668        ideal gb;
1669        gb=kStd(inputIdeal,NULL,testHomog,NULL);
1670        idNorm(gb);
1671        idSkipZeroes(gb);
1672        this->gcBasis=gb;       //write the GB into gcBasis
1673        test=save;
1674}//void getGB
1675               
1676/** \brief Compute the negative of a given intvec
1677        */             
1678inline intvec *gcone::ivNeg(const intvec *iv)
1679{
1680        //NOTE: Can't this be done without new?
1681        intvec *res;// = new intvec(iv->length());
1682        res=ivCopy(iv);
1683        *res *= (int)-1;                                               
1684        return res;
1685}
1686
1687
1688/** \brief Compute the dot product of two intvecs
1689*
1690*/
1691inline int gcone::dotProduct(const intvec &iva, const intvec &ivb)                             
1692{                       
1693        int res=0;
1694        for (int i=0;i<this->numVars;i++)
1695        {
1696                res = res+(iva[i]*ivb[i]);
1697        }
1698        return res;
1699}//int dotProduct
1700
1701/** \brief Check whether two intvecs are parallel
1702 *
1703 * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1704 */
1705inline bool gcone::isParallel(const intvec &a, const intvec &b)
1706{                       
1707        int lhs,rhs;
1708        bool res;
1709        lhs=dotProduct(a,b)*dotProduct(a,b);
1710        rhs=dotProduct(a,a)*dotProduct(b,b);
1711                        //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1712        if (lhs==rhs)
1713        {
1714                res = TRUE;
1715        }
1716        else
1717        {
1718                res = FALSE;
1719        }
1720        return res;
1721}//bool isParallel
1722               
1723/** \brief Compute an interior point of a given cone
1724 * Result will be written into intvec iv.
1725 * Any rational point is automatically converted into an integer.
1726 */
1727inline void gcone::interiorPoint( dd_MatrixPtr &M, intvec &iv) //no const &M here since we want to remove redundant rows
1728{
1729        dd_LPPtr lp,lpInt;
1730        dd_ErrorType err=dd_NoError;
1731        dd_LPSolverType solver=dd_DualSimplex;
1732        dd_LPSolutionPtr lpSol=NULL;
1733        dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
1734        dd_rowindex ddnewpos;
1735        dd_NumberType numb;     
1736        //M->representation=dd_Inequality;
1737                       
1738        //NOTE: Make this n-dimensional!
1739        //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1740                                                                       
1741        /*NOTE: Leave the following line commented out!
1742        * Otherwise it will slow down computations a great deal
1743        * */
1744//      dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
1745        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1746        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1747        int jj=1;
1748        for (int ii=0;ii<this->numVars;ii++)
1749        {
1750                dd_set_si(posRestr->matrix[ii][jj],1);
1751                jj++;                                                   
1752        }
1753        dd_MatrixAppendTo(&M,posRestr);
1754        dd_FreeMatrix(posRestr);
1755//      dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
1756        lp=dd_Matrix2LP(M, &err);
1757        if (err!=dd_NoError){WerrorS("Error during dd_Matrix2LP in gcone::interiorPoint");}
1758        if (lp==NULL){WerrorS("LP is NULL");}
1759#ifdef gfan_DEBUG
1760//                      dd_WriteLP(stdout,lp);
1761#endif 
1762                                       
1763        lpInt=dd_MakeLPforInteriorFinding(lp);
1764        if (err!=dd_NoError){WerrorS("Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint");}
1765#ifdef gfan_DEBUG
1766//                      dd_WriteLP(stdout,lpInt);
1767#endif                 
1768
1769//      dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1770        if (err!=dd_NoError)
1771        {
1772                WerrorS("Error during dd_FindRelativeInterior in gcone::interiorPoint");
1773                dd_WriteErrorMessages(stdout, err);
1774        }
1775                       
1776        dd_LPSolve(lpInt,solver,&err);  //This will not result in a point from the relative interior
1777//      if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}
1778                                                                       
1779        lpSol=dd_CopyLPSolution(lpInt);
1780//      if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
1781#ifdef gfan_DEBUG
1782        cout << "Interior point: ";
1783        for (int ii=1; ii<(lpSol->d)-1;ii++)
1784        {
1785                dd_WriteNumber(stdout,lpSol->sol[ii]);
1786        }
1787        cout << endl;
1788#endif
1789                       
1790        //NOTE The following strongly resembles parts of makeInt.
1791        //Maybe merge sometimes
1792        mpz_t kgV; mpz_init(kgV);
1793        mpz_set_str(kgV,"1",10);
1794        mpz_t den; mpz_init(den);
1795        mpz_t tmp; mpz_init(tmp);
1796        mpq_get_den(tmp,lpSol->sol[1]);
1797        for (int ii=1;ii<(lpSol->d)-1;ii++)
1798        {
1799                mpq_get_den(den,lpSol->sol[ii+1]);
1800                mpz_lcm(kgV,tmp,den);
1801                mpz_set(tmp, kgV);
1802        }
1803        mpq_t qkgV;
1804        mpq_init(qkgV);
1805        mpq_set_z(qkgV,kgV);
1806        for (int ii=1;ii<(lpSol->d)-1;ii++)
1807        {
1808                mpq_t product;
1809                mpq_init(product);
1810                mpq_mul(product,qkgV,lpSol->sol[ii]);
1811                iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
1812                mpq_clear(product);
1813        }
1814#ifdef gfan_DEBUG
1815//                      iv.show();
1816//                      cout << endl;
1817#endif
1818        mpq_clear(qkgV);
1819        mpz_clear(tmp);
1820        mpz_clear(den);
1821        mpz_clear(kgV);                 
1822                       
1823        dd_FreeLPSolution(lpSol);
1824        dd_FreeLPData(lpInt);
1825        dd_FreeLPData(lp);
1826//      set_free(ddlinset);
1827//      set_free(ddredrows);   
1828                       
1829}//void interiorPoint(dd_MatrixPtr const &M)
1830
1831inline void gcone::interiorPoint2(const dd_MatrixPtr &M, intvec &iv)
1832{
1833        KMatrix<Rational> mMat(M->rowsize+1,M->colsize);
1834        for(int ii=0;ii<M->rowsize;ii++)
1835        {
1836                for(int jj=0;jj<M->colsize-1;jj++)
1837                {
1838                        if(mpq_sgn(M->matrix[ii][jj+1])<-1)
1839                        {
1840                                mMat.set(ii,jj,-(Rational)mpq_get_d(M->matrix[ii][jj+1]));
1841                        }
1842                        else
1843                                mMat.set(ii,jj,(Rational)mpq_get_d(M->matrix[ii][jj+1]));
1844                               
1845//                      mMat.set(ii,jj,&(M->matrix[ii][jj+1]) );
1846                        cout << mpq_get_d(M->matrix[ii][jj+1]) << " ";
1847//                      int val=(int)mMat.get(ii,jj);
1848//                      cout << ii << "," << jj << endl;;
1849//                      mpq_out_str (NULL, 10, (__mpq_struct)mMat.get(ii,jj));
1850                }
1851                cout << endl;
1852                mMat.set(ii,M->colsize-1,1);
1853        }
1854        dd_WriteMatrix(stdout,M);
1855//      for(int ii=0;ii<M->rowsize;ii++)
1856//      {
1857//              cout << mMat.get(ii,ii+M->colsize) << " ";
1858//              if((ii+M->colsize)%M->colsize==0)
1859//                      cout << endl;
1860//      }
1861       
1862        Rational* mSol;
1863        int rank;
1864        int c;
1865//      dd_WriteMatrix(stdout,M);
1866        rank=mMat.solve(&mSol,&c);
1867//      for(int ii=0;ii<c;ii++)
1868//              iv[ii]=mSol[ii];
1869//              cout << mSol[ii].get_den() << "/" << mSol[ii].get_num() << endl;
1870        int gcd=1;
1871        for(int ii=0;ii<c-1;ii++)
1872                gcd += intgcd(mSol[ii].get_den_si(),mSol[ii+1].get_den_si());
1873        cout << gcd << endl;
1874        for(int ii=0;ii<iv.length();ii++)
1875                iv[ii]=(int)((mSol[ii].get_num_si()*gcd)/mSol[ii].get_den_si());
1876
1877}
1878       
1879       
1880/** \brief Copy a ring and add a weighted ordering in first place
1881 *
1882 */
1883ring gcone::rCopyAndAddWeight(const ring &r, const intvec *ivw)                         
1884{
1885        ring res=rCopy0(r);
1886        int jj;
1887                       
1888        omFree(res->order);
1889        res->order =(int *)omAlloc0(4*sizeof(int));
1890        omFree(res->block0);
1891        res->block0=(int *)omAlloc0(4*sizeof(int));
1892        omFree(res->block1);
1893        res->block1=(int *)omAlloc0(4*sizeof(int));
1894        omfree(res->wvhdl);
1895        res->wvhdl =(int **)omAlloc0(4*sizeof(int**));
1896                       
1897        res->order[0]=ringorder_a;
1898        res->block0[0]=1;
1899        res->block1[0]=res->N;;
1900        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                   
1901        res->block0[1]=1;
1902        res->block1[1]=res->N;;
1903        res->order[2]=ringorder_C;
1904
1905        int length=ivw->length();
1906        int *A=(int *)omAlloc0(length*sizeof(int));
1907        for (jj=0;jj<length;jj++)
1908        {                               
1909                A[jj]=(*ivw)[jj];                               
1910        }                       
1911        res->wvhdl[0]=(int *)A;
1912        res->block1[0]=length;
1913                       
1914        rComplete(res);
1915        return res;
1916}//rCopyAndAdd
1917               
1918ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1919{
1920        ring res=rCopy0(currRing);
1921        rComplete(res);
1922        rSetWeightVec(res,(int64*)ivw);
1923        //rChangeCurrRing(rnew);
1924        return res;
1925}
1926               
1927/** \brief Checks whether a given facet is a search facet
1928 * Determines whether a given facet of a cone is the search facet of a neighbouring cone
1929 * This is done in the following way:
1930 * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1931 * that is first crossed during the generic walk.
1932 * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1933 * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1934*/
1935//removed with r12286
1936               
1937/** \brief Check for equality of two intvecs
1938 */
1939inline bool gcone::areEqual(const intvec &a, const intvec &b)
1940{
1941        bool res=TRUE;
1942        for(int ii=0;ii<this->numVars;ii++)
1943        {
1944                if(a[ii]!=b[ii])
1945                {
1946                        res=FALSE;
1947                        break;
1948                }
1949        }
1950        return res;
1951}
1952               
1953/** \brief The reverse search algorithm
1954 */
1955//removed with r12286
1956               
1957/** \brief The new method of Markwig and Jensen
1958 * Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1959 * In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1960 * Keep a list of facets as a linked list containing an intvec and an integer matrix.
1961 * Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1962 * soon as we compute this facet again. Comparison of facets is done by...
1963 */
1964void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
1965{       
1966        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
1967        facet *SearchListAct;
1968        SearchListAct = NULL;
1969        //SearchListAct = SearchListRoot;
1970                       
1971        gcone *gcAct;
1972        gcAct = &gcRoot;
1973        gcone *gcPtr;   //Pointer to end of linked list of cones
1974        gcPtr = &gcRoot;
1975        gcone *gcNext;  //Pointer to next cone to be visited
1976        gcNext = &gcRoot;
1977        gcone *gcHead;
1978        gcHead = &gcRoot;
1979                       
1980        facet *fAct;
1981        fAct = gcAct->facetPtr;                 
1982                       
1983        ring rAct;
1984        rAct = currRing;
1985                                               
1986        int UCNcounter=gcAct->getUCN();
1987       
1988        /* Use pure SLA or writeCone2File */
1989//      enum heuristic {
1990//              ram,
1991//              hdd
1992//      };
1993//      heuristic h;
1994//      h=hdd;
1995       
1996#ifdef gfan_DEBUG
1997        cout << "NoRevs" << endl;
1998        cout << "Facets are:" << endl;
1999        gcAct->showFacets();
2000#endif                 
2001                       
2002        gcAct->getCodim2Normals(*gcAct);
2003                               
2004        //Compute unique representation of codim-2-facets
2005        gcAct->normalize();
2006                       
2007        /* Make a copy of the facet list of first cone
2008           Since the operations getCodim2Normals and normalize affect the facets
2009           we must not memcpy them before these ops!
2010        */
2011       
2012        facet *codim2Act; codim2Act = NULL;                     
2013        facet *sl2Root; //sl2Root = new facet(2);
2014        facet *sl2Act;  //sl2Act = sl2Root;
2015                       
2016        for(int ii=0;ii<this->numFacets;ii++)
2017        {
2018                //only copy flippable facets! NOTE: We do not need flipRing or any such information.
2019                if(fAct->isFlippable==TRUE)
2020                {
2021                        intvec *fNormal;
2022                        fNormal = fAct->getFacetNormal();
2023                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
2024                        {                                               
2025                                SearchListAct = SearchListRoot;                         
2026                        }
2027                        else
2028                        {                       
2029                                SearchListAct->next = new facet();
2030                                SearchListAct = SearchListAct->next;                                           
2031                        }
2032                        SearchListAct->setFacetNormal(fNormal);
2033                        SearchListAct->setUCN(this->getUCN());
2034                        SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
2035                        SearchListAct->isFlippable=TRUE;
2036                        //Copy codim2-facets                           
2037                        codim2Act=fAct->codim2Ptr;
2038                        SearchListAct->codim2Ptr = new facet(2);
2039                        sl2Root = SearchListAct->codim2Ptr;
2040                        sl2Act = sl2Root;
2041                                        //while(codim2Act!=NULL)
2042                        for(int jj=0;jj<fAct->numCodim2Facets;jj++)
2043                        {
2044                                intvec *f2Normal;
2045                                f2Normal = codim2Act->getFacetNormal();
2046                                if(jj==0)
2047                                {                                               
2048                                        sl2Act = sl2Root;
2049                                        sl2Act->setFacetNormal(f2Normal);
2050                                }
2051                                else
2052                                {
2053                                        sl2Act->next = new facet(2);
2054                                        sl2Act = sl2Act->next;
2055                                        sl2Act->setFacetNormal(f2Normal);
2056                                }
2057                                delete f2Normal;                               
2058                                codim2Act = codim2Act->next;
2059                        }
2060                        fAct = fAct->next;
2061                        delete fNormal;                 
2062                }//if(fAct->isFlippable==TRUE)                 
2063                else {fAct = fAct->next;}
2064        }//End of copying facets into SLA
2065       
2066        SearchListAct = SearchListRoot; //Set to beginning of list
2067        /*Make SearchList doubly linked*/
2068        while(SearchListAct!=NULL)
2069        {
2070                if(SearchListAct->next!=NULL)
2071                {
2072                        SearchListAct->next->prev = SearchListAct;                                     
2073                }
2074                SearchListAct = SearchListAct->next;
2075        }
2076        SearchListAct = SearchListRoot; //Set to beginning of List
2077       
2078        fAct = gcAct->facetPtr;
2079        //NOTE Disabled until code works as expected. Reenabled 2.11.2009
2080        gcAct->writeConeToFile(*gcAct);                 
2081        /*End of initialisation*/
2082       
2083        fAct = SearchListAct;
2084        /*2nd step
2085          Choose a facet from fListPtr, flip it and forget the previous cone
2086          We always choose the first facet from fListPtr as facet to be flipped
2087        */                     
2088        while((SearchListAct!=NULL))// && counter<10)
2089        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
2090                fAct = SearchListAct;
2091                               
2092                while(fAct!=NULL)
2093//              while( (fAct->getUCN() == fAct->next->getUCN()) )               
2094                {       //Since SLA should only contain flippables there should be no need to check for that                   
2095                        gcAct->flip(gcAct->gcBasis,fAct);                       
2096                        //NOTE rCopy needed?
2097                        ring rTmp=rCopy(fAct->flipRing);
2098//                      ring rTmp=fAct->flipRing; //segfaults
2099                        rComplete(rTmp);                       
2100                        rChangeCurrRing(rTmp);
2101                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);//copy constructor!
2102                        /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
2103                         * Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
2104                         * NOTE: Can this cause trouble with the destructor?
2105                         * Answer: Yes it bloody well does!
2106                         * However I'd like to delete this data here, since if we have a cone with many many facets it
2107                         * should pay to get rid of the flipGB as soon as possible.
2108                         * Destructor must be equipped with necessary checks.
2109                        */
2110                        idDelete((ideal *)&fAct->flipGB);
2111                        rDelete(fAct->flipRing);
2112                       
2113                        gcTmp->getConeNormals(gcTmp->gcBasis, FALSE);   
2114                        gcTmp->getCodim2Normals(*gcTmp);
2115                        gcTmp->normalize();
2116                        //gcTmp->ddFacets should not be needed anymore, so
2117                        dd_FreeMatrix(gcTmp->ddFacets);
2118#ifdef gfan_DEBUG
2119//                      gcTmp->showFacets(1);
2120#endif
2121                        /*add facets to SLA here*/
2122                        SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
2123                        if(gfanHeuristic==1)
2124                        {
2125                                gcTmp->writeConeToFile(*gcTmp);
2126                                //The for loop is no longer needed
2127//                              for(int ii=0;ii<IDELEMS(gcTmp->gcBasis);ii++)
2128//                              {
2129//                                      pDelete(&gcTmp->gcBasis->m[ii]);
2130//                              }
2131                                idDelete((ideal*)&gcTmp->gcBasis);//Whonder why?
2132                                //If you use the following make sure it is uncommented in readConeFromFile
2133//                              rDelete(gcTmp->baseRing);
2134                        }                       
2135#ifdef gfan_DEBUG
2136//                      if(SearchListRoot!=NULL)
2137//                              gcTmp->showSLA(*SearchListRoot);
2138#endif
2139                        rChangeCurrRing(gcAct->baseRing);
2140                        rDelete(rTmp);
2141                        //doubly linked for easier removal
2142                        gcTmp->prev = gcPtr;
2143                        gcPtr->next=gcTmp;
2144                        gcPtr=gcPtr->next;
2145                        if(fAct->getUCN() == fAct->next->getUCN())
2146                        {
2147                                fAct=fAct->next;
2148                        }
2149                        else
2150                                break;
2151//                      fAct=fAct->next;
2152                }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );         
2153                //Search for cone with smallest UCN
2154                gcNext = gcHead;
2155                while(gcNext!=NULL && SearchListRoot!=NULL && gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && gcNext!=(gcone * const)0xfbfbfbfb)
2156                {                       
2157                        if( gcNext->getUCN() == SearchListRoot->getUCN() )
2158                        {//NOTE: Implement heuristic to be used!
2159                                if(gfanHeuristic==0)
2160                                {
2161                                        gcAct = gcNext;
2162                                        //Seems better to not to use rCopy here
2163//                                      rAct=rCopy(gcAct->baseRing);
2164                                        rAct=gcAct->baseRing;
2165                                        rComplete(rAct);
2166                                        rChangeCurrRing(rAct);                                         
2167                                        break;
2168                                }
2169                                else if(gfanHeuristic==1)
2170                                {
2171                                        gcone *gcDel;
2172                                        gcDel = gcAct;                                 
2173                                        gcAct = gcNext;
2174//                                      rDelete(gcDel->baseRing);
2175                                        //Read st00f from file
2176                                        //implant the GB into gcAct
2177                                        readConeFromFile(gcAct->getUCN(), gcAct);
2178//                                      rAct=rCopy(gcAct->baseRing);
2179                                        rAct=gcAct->baseRing;
2180                                        rComplete(rAct);
2181                                        rChangeCurrRing(rAct);
2182                                        break;
2183                                }                               
2184                        }                       
2185                        gcNext = gcNext->next;
2186                }
2187                UCNcounter++;
2188                //SearchListAct = SearchListAct->next;
2189                //SearchListAct = fAct->next;
2190                SearchListAct = SearchListRoot;
2191        }
2192        cout << endl << "Found " << counter << " cones - terminating" << endl;
2193}//void noRevS(gcone &gc)       
2194               
2195               
2196/** \brief Make a set of rational vectors into integers
2197 *
2198 * We compute the lcm of the denominators and multiply with this to get integer values.         
2199 * \param dd_MatrixPtr,intvec
2200 */
2201inline void gcone::makeInt(const dd_MatrixPtr &M, const int line, intvec &n)
2202{                       
2203//      mpz_t denom[this->numVars];
2204        mpz_t *denom = new mpz_t[this->numVars];
2205        for(int ii=0;ii<this->numVars;ii++)
2206        {
2207                mpz_init_set_str(denom[ii],"0",10);
2208        }
2209               
2210        mpz_t kgV,tmp;
2211        mpz_init(kgV);
2212        mpz_init(tmp);
2213                       
2214        for (int ii=0;ii<(M->colsize)-1;ii++)
2215        {
2216                mpz_t z;
2217                mpz_init(z);
2218                mpq_get_den(z,M->matrix[line-1][ii+1]);
2219                mpz_set( denom[ii], z);
2220                mpz_clear(z);                           
2221        }
2222                                               
2223        /*Compute lcm of the denominators*/
2224        mpz_set(tmp,denom[0]);
2225        for (int ii=0;ii<(M->colsize)-1;ii++)
2226        {
2227                mpz_lcm(kgV,tmp,denom[ii]);
2228                mpz_set(tmp,kgV);                               
2229        }
2230        mpz_clear(tmp); 
2231        /*Multiply the nominators by kgV*/
2232        mpq_t qkgV,res;
2233        mpq_init(qkgV);
2234        mpq_set_str(qkgV,"1",10);
2235        mpq_canonicalize(qkgV);
2236                       
2237        mpq_init(res);
2238        mpq_set_str(res,"1",10);
2239        mpq_canonicalize(res);
2240                       
2241        mpq_set_num(qkgV,kgV);
2242                       
2243//                      mpq_canonicalize(qkgV);
2244        for (int ii=0;ii<(M->colsize)-1;ii++)
2245        {
2246                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
2247                n[ii]=(int)mpz_get_d(mpq_numref(res));
2248        }
2249        delete [] denom;
2250        mpz_clear(kgV);
2251        mpq_clear(qkgV); mpq_clear(res);                       
2252                       
2253}
2254/**
2255 * For all codim-2-facets we compute the gcd of the components of the facet normal and
2256 * divide it out. Thus we get a normalized representation of each
2257 * (codim-2)-facet normal, i.e. a primitive vector
2258 */
2259inline void gcone::normalize()
2260{
2261        int *ggT = new int;
2262                *ggT=1;
2263        facet *fAct;
2264        facet *codim2Act;
2265        fAct = this->facetPtr;
2266        codim2Act = fAct->codim2Ptr;
2267//      intvec *n = new intvec(this->numVars);
2268                       
2269        //while(codim2Act->next!=NULL)
2270        while(fAct!=NULL)
2271        {
2272                while(codim2Act!=NULL)
2273                {                               
2274                        intvec *n;
2275                        n=codim2Act->getFacetNormal();
2276                        for(int ii=0;ii<this->numVars;ii++)
2277                        {
2278                                *ggT = intgcd((*ggT),(*n)[ii]);
2279                        }
2280                        for(int ii=0;ii<this->numVars;ii++)
2281                        {
2282                                (*n)[ii] = ((*n)[ii])/(*ggT);
2283                        }
2284                        codim2Act->setFacetNormal(n);
2285                        codim2Act = codim2Act->next;
2286                        delete n;
2287                }
2288                fAct = fAct->next;
2289        }
2290        delete ggT;
2291//      delete n;
2292                               
2293}
2294/** \brief Enqueue new facets into the searchlist
2295 * The searchlist (SLA for short) is implemented as a FIFO
2296 * Checks are done as follows:
2297 * 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
2298 * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
2299 *      facet from SLA and do not add fAct.
2300 * It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
2301 * can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA.
2302 * Takes ptr to search list root
2303 * Returns a pointer to new first element of Searchlist
2304 */
2305facet * gcone::enqueueNewFacets(facet *f)
2306{
2307#ifdef gfanp
2308        timeval start, end;
2309        gettimeofday(&start, 0);
2310#endif
2311        facet *slHead;
2312        slHead = f;
2313        facet *slAct;   //called with f=SearchListRoot
2314        slAct = f;
2315        facet *slEnd;   //Pointer to end of SLA
2316        slEnd = f;
2317//      facet *slEndStatic;     //marks the end before a new facet is added             
2318        facet *fAct;
2319        fAct = this->facetPtr;
2320        facet *codim2Act;
2321        codim2Act = this->facetPtr->codim2Ptr;
2322        facet *sl2Act;
2323        sl2Act = f->codim2Ptr;
2324        /** Pointer to a facet that will be deleted */
2325        volatile facet *deleteMarker;
2326        deleteMarker = NULL;
2327                       
2328        /** Flag to indicate a facet that should be added to SLA*/
2329//      bool doNotAdd=FALSE;
2330        /** \brief  Flag to mark a facet that might be added
2331         * The following case may occur:
2332         * A facet fAct is found to be parallel but not equal to the current facet slAct from SLA.
2333         * This does however not mean that there does not exist a facet behind the current slAct that is actually equal
2334         * to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
2335         * FALSE and the according slAct is deleted.
2336         * If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
2337         */
2338//      volatile bool maybe=FALSE;
2339        /**A facet was removed, lengthOfSearchlist-- thus we must not rely on
2340         * if(notParallelCtr==lengthOfSearchList) but rather
2341         * if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
2342         */
2343        volatile bool removalOccured=FALSE;
2344        int ctr=0;      //encountered equalities in SLA
2345        int notParallelCtr=0;
2346        int lengthOfSearchList=1;
2347        while(slEnd->next!=NULL)
2348        {
2349                slEnd=slEnd->next;
2350                lengthOfSearchList++;
2351        }
2352        /*1st step: compare facetNormals*/
2353//      intvec *fNormal=NULL;
2354//      intvec *slNormal=NULL;
2355                       
2356        while(fAct!=NULL)
2357        {                                               
2358                if(fAct->isFlippable==TRUE)
2359                {
2360                        intvec *fNormal=NULL;
2361//                      intvec *slNormal=NULL;
2362                        fNormal=fAct->getFacetNormal();
2363                        slAct = slHead;
2364                        notParallelCtr=0;
2365                        /*If slAct==NULL and fAct!=NULL
2366                        we just copy all remaining facets into SLA*/
2367                        if(slAct==NULL)
2368                        {
2369                                facet *fCopy;
2370                                fCopy = fAct;
2371                                while(fCopy!=NULL)
2372                                {
2373                                        if(slAct==NULL)
2374                                        {
2375                                                slAct = new facet(*fCopy);
2376                                                slHead = slAct;                                                         
2377                                        }
2378                                        else
2379                                        {
2380                                                slAct->next = new facet(*fCopy);
2381                                                slAct = slAct->next;
2382                                        }
2383                                        fCopy = fCopy->next;
2384                                }
2385                                break;
2386                        }
2387                        /*End of dumping into SLA*/
2388                        while(slAct!=NULL)
2389                        //while(slAct!=slEndStatic->next)
2390                        {
2391                                intvec *slNormal=NULL;
2392                                removalOccured=FALSE;
2393                                slNormal = slAct->getFacetNormal();
2394#ifdef gfan_DEBUG
2395                                cout << "Checking facet (";
2396                                fNormal->show(1,1);
2397                                cout << ") against (";
2398                                slNormal->show(1,1);
2399                                cout << ")" << endl;
2400#endif                         
2401                                if(areEqual(fAct,slAct))
2402                                {
2403                                        deleteMarker = slAct;
2404                                        if(slAct==slHead)
2405                                        {                                               
2406                                                slHead = slAct->next;                                           
2407                                                if(slHead!=NULL)
2408                                                        slHead->prev = NULL;                                           
2409                                        }
2410                                        else if (slAct==slEnd)
2411                                        {
2412                                                slEnd=slEnd->prev;
2413                                                slEnd->next = NULL;                                             
2414                                        }                                                               
2415                                        else
2416                                        {
2417                                                slAct->prev->next = slAct->next;
2418                                                slAct->next->prev = slAct->prev;                                               
2419                                        }
2420                                        removalOccured=TRUE;
2421                                        lengthOfSearchList--;
2422                                        if(deleteMarker!=NULL)
2423                                        {
2424//                                              delete deleteMarker;
2425//                                              deleteMarker=NULL;
2426                                        }
2427#ifdef gfan_DEBUG
2428                                        cout << "Removing (";
2429                                        fNormal->show(1,1);
2430                                        cout << ") from list" << endl;
2431#endif
2432                                        delete slNormal;
2433                                        break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
2434                                }               
2435                               
2436                                slAct = slAct->next;
2437                                /* NOTE The following lines must not be here but rather called inside the if-block above.
2438                                Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now,
2439                                (not nice!) but since they are in seperate branches of the if-statement there should not
2440                                be a way it gets called twice thus ommiting one facet:
2441                                slAct = slAct->next;*/
2442                                if(deleteMarker!=NULL)
2443                                {                                               
2444//                                      delete deleteMarker;
2445//                                      deleteMarker=NULL;
2446                                }
2447                                delete slNormal;
2448                                                //if slAct was marked as to be deleted, delete it here!
2449                        }//while(slAct!=NULL)
2450                        if(removalOccured==FALSE)
2451//                      if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) || (doNotAdd==FALSE) )
2452                                        //if( (notParallelCtr==lengthOfSearchList ) || doNotAdd==FALSE )
2453                        {
2454#ifdef gfan_DEBUG
2455                                cout << "Adding facet (";
2456                                fNormal->show(1,0);
2457                                cout << ") to SLA " << endl;
2458#endif
2459                                                //Add fAct to SLA
2460                                facet *marker;
2461                                marker = slEnd;
2462                                facet *f2Act;
2463                                f2Act = fAct->codim2Ptr;
2464                                               
2465                                slEnd->next = new facet();
2466                                slEnd = slEnd->next;
2467                                facet *slEndCodim2Root;
2468                                facet *slEndCodim2Act;
2469                                slEndCodim2Root = slEnd->codim2Ptr;
2470                                slEndCodim2Act = slEndCodim2Root;
2471                                                               
2472                                slEnd->setUCN(this->getUCN());
2473                                slEnd->isFlippable = TRUE;
2474                                slEnd->setFacetNormal(fNormal);
2475                                slEnd->prev = marker;
2476                                //Copy codim2-facets
2477//                              intvec *f2Normal=new intvec(this->numVars);
2478                                while(f2Act!=NULL)
2479                                {
2480                                        intvec *f2Normal;
2481                                        f2Normal=f2Act->getFacetNormal();
2482                                        if(slEndCodim2Root==NULL)
2483                                        {
2484                                                slEndCodim2Root = new facet(2);
2485                                                slEnd->codim2Ptr = slEndCodim2Root;                     
2486                                                slEndCodim2Root->setFacetNormal(f2Normal);
2487                                                slEndCodim2Act = slEndCodim2Root;
2488                                        }
2489                                        else                                   
2490                                        {
2491                                                slEndCodim2Act->next = new facet(2);
2492                                                slEndCodim2Act = slEndCodim2Act->next;
2493                                                slEndCodim2Act->setFacetNormal(f2Normal);
2494                                        }
2495                                        f2Act = f2Act->next;
2496                                        delete f2Normal;
2497                                }
2498                                lengthOfSearchList++;                           
2499                        }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
2500                        fAct = fAct->next;
2501                        delete fNormal;
2502//                      delete slNormal;
2503                }//if(fAct->isFlippable==TRUE)
2504                else
2505                {
2506                        fAct = fAct->next;
2507                }
2508        }//while(fAct!=NULL)
2509#ifdef gfanp
2510        gettimeofday(&end, 0);
2511        time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2512#endif                                         
2513        return slHead;
2514}//addC2N
2515               
2516/** \brief Compute the gcd of two ints
2517 */
2518inline int gcone::intgcd(const int a, const int b)
2519{
2520        int r, p=a, q=b;
2521        if(p < 0)
2522                p = -p;
2523        if(q < 0)
2524                q = -q;
2525        while(q != 0)
2526        {
2527                r = p % q;
2528                p = q;
2529                q = r;
2530        }
2531        return p;
2532}               
2533               
2534/** \brief Construct a dd_MatrixPtr from a cone's list of facets
2535 * NO LONGER USED
2536 */
2537inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
2538{
2539        facet *fAct;
2540        fAct = gc.facetPtr;
2541       
2542        dd_MatrixPtr M;
2543        dd_rowrange ddrows;
2544        dd_colrange ddcols;
2545        ddcols=(this->numVars)+1;
2546        ddrows=this->numFacets;
2547        dd_NumberType numb = dd_Integer;
2548        M=dd_CreateMatrix(ddrows,ddcols);                       
2549                       
2550        int jj=0;
2551       
2552        while(fAct!=NULL)
2553        {
2554                intvec *comp;
2555                comp = fAct->getFacetNormal();
2556                for(int ii=0;ii<this->numVars;ii++)
2557                {
2558                        dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
2559                }
2560                jj++;
2561                delete comp;
2562                fAct=fAct->next;                               
2563        }                       
2564        return M;
2565}
2566               
2567/** \brief Write information about a cone into a file on disk
2568 *
2569 * This methods writes the information needed for the "second" method into a file.
2570 * The file's is divided in sections containing information on
2571 * 1) the ring
2572 * 2) the cone's Groebner Basis
2573 * 3) the cone's facets
2574 * Each line contains exactly one date
2575 * Each section starts with its name in CAPITALS
2576 */
2577inline void gcone::writeConeToFile(const gcone &gc, bool usingIntPoints)
2578{
2579        int UCN=gc.UCN;
2580        stringstream ss;
2581        ss << UCN;
2582        string UCNstr = ss.str();               
2583                       
2584        string prefix="/tmp/cone";
2585//      string prefix="./";     //crude hack -> run tests in separate directories
2586        string suffix=".sg";
2587        string filename;
2588        filename.append(prefix);
2589        filename.append(UCNstr);
2590        filename.append(suffix);
2591       
2592       
2593//      int thisPID = getpid();
2594//      ss << thisPID;
2595//      string strPID = ss.str();
2596//      string prefix="./";
2597                                               
2598        ofstream gcOutputFile(filename.c_str());
2599        facet *fAct;
2600        fAct = gc.facetPtr;                     
2601        facet *f2Act;
2602        f2Act=fAct->codim2Ptr;
2603       
2604        char *ringString = rString(currRing);
2605                       
2606        if (!gcOutputFile)
2607        {
2608                cout << "Error opening file for writing in writeConeToFile" << endl;
2609        }
2610        else
2611        {       
2612                gcOutputFile << "UCN" << endl;
2613                gcOutputFile << this->UCN << endl;
2614                gcOutputFile << "RING" << endl; 
2615                gcOutputFile << ringString << endl;
2616                gcOutputFile << "GCBASISLENGTH" << endl;
2617                gcOutputFile << IDELEMS(gc.gcBasis) << endl;                   
2618                //Write this->gcBasis into file
2619                gcOutputFile << "GCBASIS" << endl;                             
2620                for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
2621                {                                       
2622                        gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
2623                }                               
2624                                       
2625                gcOutputFile << "FACETS" << endl;                                                               
2626               
2627                while(fAct!=NULL)
2628                {       
2629                        intvec *iv;
2630                        iv=fAct->getFacetNormal();
2631                        f2Act=fAct->codim2Ptr;
2632                        for (int ii=0;ii<iv->length();ii++)
2633                        {
2634                                if (ii<iv->length()-1)
2635                                {
2636                                        gcOutputFile << (*iv)[ii] << ",";
2637                                }
2638                                else
2639                                {
2640                                        gcOutputFile << (*iv)[ii] << " ";
2641                                }
2642                        }
2643                        delete iv;
2644                        while(f2Act!=NULL)
2645                        {
2646                                intvec *iv2;
2647                                iv2=f2Act->getFacetNormal();   
2648                                for(int jj=0;jj<iv2->length();jj++)
2649                                {
2650                                        if (jj<iv2->length()-1)
2651                                        {
2652                                                gcOutputFile << (*iv2)[jj] << ",";
2653                                        }
2654                                        else
2655                                        {
2656                                                gcOutputFile << (*iv2)[jj] << " ";
2657                                        }
2658                                }
2659                                delete iv2;
2660                                f2Act = f2Act->next;
2661                        }
2662                        gcOutputFile << endl;
2663                        fAct=fAct->next;
2664                }                       
2665        gcOutputFile.close();
2666        }
2667                       
2668}//writeConeToFile(gcone const &gc)
2669               
2670/** \brief Reads a cone from a file identified by its number
2671 */
2672inline void gcone::readConeFromFile(int UCN, gcone *gc)
2673{
2674        //int UCN=gc.UCN;
2675        stringstream ss;
2676        ss << UCN;
2677        string UCNstr = ss.str();
2678//      string line;
2679//      string strGcBasisLength;
2680//      string strMonom, strCoeff, strCoeffNom, strCoeffDenom;
2681        int gcBasisLength=0;
2682//      int intCoeff=1;
2683//      int intCoeffNom=1;              //better (number) but that's not compatible with stringstream;
2684//      int intCoeffDenom=1;
2685       
2686//      bool hasCoeffInQ = FALSE;       //does the polynomial have rational coeff?
2687//      bool hasNegCoeff = FALSE;       //or a negative one?
2688        size_t found;                   //used for find_first_(not)_of
2689
2690        string prefix="/tmp/cone";
2691        string suffix=".sg";
2692        string filename;
2693        filename.append(prefix);
2694        filename.append(UCNstr);
2695        filename.append(suffix);
2696                                       
2697        ifstream gcInputFile(filename.c_str(), ifstream::in);
2698       
2699        ring saveRing=currRing; 
2700        //Comment the following if you uncomment the if(line=="RING") part below
2701        rChangeCurrRing(gc->baseRing);
2702       
2703        while( !gcInputFile.eof() )
2704        {
2705                string line;
2706                getline(gcInputFile,line);
2707//              hasCoeffInQ = FALSE;
2708//              hasNegCoeff = FALSE;
2709               
2710                if(line=="RING")
2711                {
2712//                      getline(gcInputFile,line);
2713//                      found = line.find("a(");
2714//                      line.erase(0,found+2);
2715//                      string strweight;
2716//                      strweight=line.substr(0,line.find_first_of(")"));
2717//                      intvec *iv=new intvec(this->numVars);
2718//                      for(int ii=0;ii<this->numVars;ii++)
2719//                      {
2720//                              string weight;
2721//                              weight=line.substr(0,line.find_first_of(",)"));                         
2722//                              (*iv)[ii]=atoi(weight.c_str());
2723//                              line.erase(0,line.find_first_of(",)")+1);
2724//                      }
2725//                      ring newRing;
2726//                      if(currRing->order[0]!=ringorder_a)
2727//                      {
2728//                              newRing=rCopyAndAddWeight(currRing,iv);
2729//                      }
2730//                      else
2731//                      {                       
2732//                              newRing=rCopy0(currRing);
2733//                              int length=this->numVars;
2734//                              int *A=(int *)omAlloc0(length*sizeof(int));
2735//                              for(int jj=0;jj<length;jj++)
2736//                              {
2737//                                      A[jj]=-(*iv)[jj];
2738//                              }
2739//                              omFree(newRing->wvhdl[0]);
2740//                              newRing->wvhdl[0]=(int*)A;
2741//                              newRing->block1[0]=length;
2742//                      }
2743//                      delete iv;
2744//                      rComplete(newRing);
2745//                      gc->baseRing=rCopy(newRing);
2746//                      if(currRing!=gc->baseRing)
2747//                              rChangeCurrRing(gc->baseRing);
2748                }
2749               
2750                if(line=="GCBASISLENGTH")
2751                {
2752                        string strGcBasisLength;
2753                        getline(gcInputFile, line);
2754                        strGcBasisLength = line;
2755                        int size=atoi(strGcBasisLength.c_str());
2756                        gcBasisLength=size;             
2757                        gc->gcBasis=idInit(size,1);
2758                }
2759                if(line=="GCBASIS")
2760                {
2761                        for(int jj=0;jj<gcBasisLength;jj++)
2762                        {
2763                                getline(gcInputFile,line);
2764                                //magically convert strings into polynomials
2765                                //polys.cc:p_Read
2766                                //check until first occurance of + or -
2767                                //data or c_str
2768                                poly strPoly=pInit();//Ought to be inside the while loop, but that will eat your memory
2769                                poly resPoly=pInit();   //The poly to be read in                                                       
2770                                while(!line.empty())
2771                                {
2772//                                      poly strPoly=pInit();
2773                                        number nCoeff=nInit(1);
2774                                        number nCoeffNom=nInit(1);
2775                                        number nCoeffDenom=nInit(1);
2776                                        string strMonom, strCoeff, strCoeffNom, strCoeffDenom;
2777                                        bool hasCoeffInQ = FALSE;       //does the polynomial have rational coeff?
2778                                        bool hasNegCoeff = FALSE;       //or a negative one?
2779                                        found = line.find_first_of("+-");       //get the first monomial
2780                                        string tmp;
2781                                        tmp=line[found];
2782//                                      if(found!=0 && (tmp.compare("-")==0) )
2783//                                              hasNegCoeff = TRUE;     //We have a coeff < 0
2784                                        if(found==0)
2785                                        {
2786                                                if(tmp.compare("-")==0)
2787                                                        hasNegCoeff = TRUE;
2788                                                line.erase(0,1);        //remove leading + or -
2789                                                found = line.find_first_of("+-");       //adjust found
2790                                        }
2791                                        strMonom = line.substr(0,found);
2792                                        line.erase(0,found);
2793                                        found = strMonom.find_first_of("/");
2794                                        if(found!=string::npos) //i.e. "/" exists in strMonom
2795                                        {
2796                                                hasCoeffInQ = TRUE;
2797                                                strCoeffNom=strMonom.substr(0,found);                                           
2798                                                strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890",found+1));
2799                                                strMonom.erase(0,found);
2800                                                strMonom.erase(0,strMonom.find_first_not_of("1234567890/"));                   
2801                                                nRead(strCoeffNom.c_str(), &nCoeffNom);
2802                                                nRead(strCoeffDenom.c_str(), &nCoeffDenom);
2803                                        }
2804                                        else
2805                                        {
2806                                                found = strMonom.find_first_not_of("1234567890");
2807                                                strCoeff = strMonom.substr(0,found);
2808                                                if(!strCoeff.empty())
2809                                                {
2810                                                        nRead(strCoeff.c_str(),&nCoeff);
2811                                                }
2812                                                else
2813                                                {
2814//                                                      intCoeff = 1;
2815                                                        nCoeff = nInit(1);
2816                                                }
2817                                                                                               
2818                                        }
2819                                        const char* monom = strMonom.c_str();
2820                                               
2821                                        p_Read(monom,strPoly,currRing); //strPoly:=monom                               
2822                                        switch (hasCoeffInQ)
2823                                        {
2824                                                case TRUE:
2825                                                        if(hasNegCoeff)
2826                                                                nCoeffNom=nNeg(nCoeffNom);
2827                                                        pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom));
2828                                                        break;
2829                                                case FALSE:
2830                                                        if(hasNegCoeff)
2831                                                                nCoeff=nNeg(nCoeff);                                                   
2832                                                        if(!nIsOne(nCoeff))
2833                                                        {
2834//                                                              if(hasNegCoeff)
2835//                                                                      intCoeff *= -1;
2836//                                                              pSetCoeff(strPoly,(number) intCoeff);
2837                                                                pSetCoeff(strPoly, nCoeff );
2838                                                        }
2839                                                        break;
2840                                                                                                       
2841                                        }
2842                                                //pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1???
2843                                        if(pIsConstantComp(resPoly))
2844                                                resPoly=pCopy(strPoly);                                                 
2845                                        else
2846                                                resPoly=pAdd(resPoly,strPoly);
2847                                        nDelete(&nCoeff);
2848                                        nDelete(&nCoeffNom);
2849                                        nDelete(&nCoeffDenom);
2850//                                      pDelete(&strPoly);
2851                                }//while(!line.empty())                 
2852                                gc->gcBasis->m[jj]=pCopy(resPoly);
2853                                pDelete(&resPoly);      //reset
2854//                              pDelete(&strPoly);      //NOTE Crashes                         
2855                        }
2856                        break;
2857                }//if(line=="GCBASIS")         
2858        }//while(!gcInputFile.eof())   
2859        gcInputFile.close();
2860        rChangeCurrRing(saveRing);
2861}
2862       
2863ring gcone::getBaseRing()
2864{
2865        return rCopy(this->baseRing);
2866}
2867/** \brief Gather the output
2868* List of lists
2869*\param *gc Pointer to gcone, preferably gcRoot ;-)
2870*\param n the number of cones
2871*
2872*/
2873lists lprepareResult(gcone *gc, int n)
2874{
2875        gcone *gcAct;
2876        gcAct = gc;     
2877        facet *fAct;
2878        fAct = gc->facetPtr;
2879       
2880        lists res=(lists)omAllocBin(slists_bin);
2881        res->Init(n);   //initialize to store n cones
2882        for(int ii=0;ii<n;ii++)
2883        {
2884                res->m[ii].rtyp=LIST_CMD;
2885                lists l=(lists)omAllocBin(slists_bin);
2886                l->Init(6);
2887                l->m[0].rtyp=INT_CMD;
2888                l->m[0].data=(void*)gcAct->getUCN();
2889                l->m[1].rtyp=IDEAL_CMD;
2890                /*The following is necessary for leaves in the tree of cones
2891                * Since we don't use them in the computation and gcBasis is
2892                * set to (poly)NULL in noRevS we need to get this back here.
2893                */
2894//              if(gcAct->gcBasis->m[0]==(poly)NULL)
2895                if(gfanHeuristic==1)
2896                        gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
2897//              ring saveRing=currRing;
2898//              ring tmpRing=gcAct->getBaseRing;
2899//              rChangeCurrRing(tmpRing);
2900//              l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
2901                l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());
2902//              rChangeCurrRing(saveRing);
2903
2904                l->m[2].rtyp=INTVEC_CMD;
2905                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
2906                l->m[2].data=(void*)ivCopy(&iv);
2907               
2908                l->m[3].rtyp=LIST_CMD;
2909                        lists lCodim2List = (lists)omAllocBin(slists_bin);
2910                        lCodim2List->Init(gcAct->numFacets); 
2911                        fAct = gcAct->facetPtr;//fAct->codim2Ptr;
2912                        int jj=0;
2913                        while(fAct!=NULL && jj<gcAct->numFacets)
2914                        {
2915                                lCodim2List->m[jj].rtyp=INTVEC_CMD;
2916                                intvec ivC2=(gcAct->f2M(gcAct,fAct,2));
2917                                lCodim2List->m[jj].data=(void*)ivCopy(&ivC2);
2918                                jj++;
2919                                fAct = fAct->next;
2920                        }               
2921                l->m[3].data=(void*)lCodim2List;               
2922                l->m[4].rtyp=RING_CMD;
2923                l->m[4].data=(void*)(gcAct->getBaseRing());             
2924                l->m[5].rtyp=INT_CMD;
2925                l->m[5].data=(void*)gcAct->getPredUCN();
2926                res->m[ii].data=(void*)l;
2927                gcAct = gcAct->next;
2928        }       
2929        return res;
2930}
2931/** \brief Write facets of a cone into a matrix
2932* Takes a pointer to a facet as 2nd arg
2933* f should always point to gc->facetPtr
2934* param n is used to determine whether it operates in codim 1 or 2
2935*
2936*/
2937inline intvec gcone::f2M(gcone *gc, facet *f, int n)
2938{
2939        facet *fAct;
2940        intvec *res;//=new intvec(this->numVars);       
2941//      int codim=n;
2942//      int bound;
2943//      if(f==gc->facetPtr)
2944        if(n==1)
2945        {
2946                intvec *m1Res=new intvec(gc->numFacets,gc->numVars,0);
2947                res = ivCopy(m1Res);
2948                fAct = gc->facetPtr;
2949                delete m1Res;
2950//              bound = gc->numFacets*(this->numVars);         
2951        }
2952        else
2953        {
2954                fAct = f->codim2Ptr;
2955                intvec *m2Res = new intvec(f->numCodim2Facets,gc->numVars,0);
2956                res = ivCopy(m2Res);
2957                delete m2Res;   
2958//              bound = fAct->numCodim2Facets*(this->numVars);
2959
2960        }
2961        int ii=0;
2962        while(fAct!=NULL )//&& ii < bound )
2963        {
2964                intvec *fNormal;
2965                fNormal = fAct->getFacetNormal();
2966                for(int jj=0;jj<this->numVars;jj++)
2967                {
2968                        (*res)[ii]=(*fNormal)[jj];
2969                        ii++;
2970                }
2971                fAct = fAct->next;
2972                delete fNormal;
2973        }       
2974        return *res;
2975}
2976
2977int gcone::counter=0;
2978int gfanHeuristic;
2979#ifdef gfanp
2980float gcone::time_getConeNormals;
2981float gcone::time_getCodim2Normals;
2982float gcone::time_flip;
2983float gcone::t_markings;
2984float gcone::t_dd;
2985float gcone::time_enqueue;
2986float gcone::time_computeInv;
2987float gcone::t_ddMC;
2988float gcone::t_mI;
2989float gcone::t_iP;
2990#endif
2991// ideal gfan(ideal inputIdeal, int h)
2992lists gfan(ideal inputIdeal, int h)
2993{
2994        lists lResList; //this is the object we return
2995       
2996        if(rHasGlobalOrdering(currRing))
2997        {       
2998        int numvar = pVariables; 
2999        gfanHeuristic = h;
3000       
3001        enum searchMethod {
3002                reverseSearch,
3003                noRevS
3004        };
3005       
3006        searchMethod method;
3007        method = noRevS;
3008
3009        ring inputRing=currRing;        // The ring the user entered
3010        ring rootRing;                  // The ring associated to the target ordering
3011        ideal res;     
3012        facet *fRoot;   
3013        dd_set_global_constants();
3014        if(method==noRevS)
3015        {
3016                gcone *gcRoot = new gcone(currRing,inputIdeal);
3017                gcone *gcAct;
3018                gcAct = gcRoot;
3019                gcAct->numVars=pVariables;
3020                gcAct->getGB(inputIdeal);               
3021#ifndef NDEBUG
3022                cout << "GB of input ideal is:" << endl;
3023//              idShow(gcAct->gcBasis);
3024#endif
3025                if(gcAct->isMonomial(gcAct->gcBasis))
3026                {
3027                        WerrorS("Monomial input - terminating");
3028                        exit(-1);
3029                        gcAct->getConeNormals(gcAct->gcBasis);
3030                        lResList=lprepareResult(gcAct,1);
3031                        dd_free_global_constants;
3032                        //This is filthy
3033                        return lResList;
3034                }
3035                cout << endl;
3036                gcAct->getConeNormals(gcAct->gcBasis);
3037                gcAct->noRevS(*gcAct);
3038                rChangeCurrRing(inputRing);
3039                //res=gcAct->gcBasis;
3040                //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
3041                res = inputIdeal; 
3042                lResList=lprepareResult(gcRoot,gcRoot->getCounter());
3043                /*Cleanup*/
3044                gcone *gcDel;   
3045                gcDel = gcRoot;
3046                gcAct = gcRoot;
3047                while(gcAct!=NULL)
3048                {
3049                        gcDel = gcAct;
3050                        gcAct = gcAct->next;
3051                        delete gcDel;
3052                }
3053        }
3054        dd_free_global_constants();
3055        }//rHasGlobalOrdering
3056        else
3057        {
3058                WerrorS("Ring has non-global ordering - terminating");
3059                lResList->Init(1);
3060                lResList->m[0].rtyp=INT_CMD;
3061                int ires=0;
3062                lResList->m[0].data=(void*)&ires;
3063        }
3064        //gcone::counter=0;
3065        /*Return result*/
3066#ifdef gfanp
3067        cout << "t_getConeNormals:" << gcone::time_getConeNormals << endl;
3068        cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl;
3069        cout << "  t_ddMC:" << gcone::t_ddMC << endl;
3070        cout << "  t_mI:" << gcone::t_mI << endl;
3071        cout << "  t_iP:" << gcone::t_iP << endl;
3072        cout << "t_Flip:" << gcone::time_flip << endl;
3073        cout << "  t_markings:" << gcone::t_markings << endl;
3074        cout << "  t_dd:" << gcone::t_dd << endl;
3075        cout << "t_computeInv:" << gcone::time_computeInv << endl;
3076        cout << "t_enqueue:" << gcone::time_enqueue << endl;
3077#endif
3078        return lResList;
3079}
3080
3081#endif
Note: See TracBrowser for help on using the repository browser.