source: git/kernel/gfan.cc @ fec53d

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