source: git/kernel/gfan.cc @ abe5c8

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