source: git/kernel/gfan.cc @ 17f35f6

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