source: git/kernel/gfan.cc @ 5ece85

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