source: git/kernel/gfan.cc @ ce41a2e

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