source: git/kernel/gfan.cc @ 706b89

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