source: git/kernel/gfan.cc @ d4f1b95

spielwiese
Last change on this file since d4f1b95 was d4f1b95, checked in by Martin Monerjan, 14 years ago
Copy constructor for facets pDelete git-svn-id: file:///usr/local/Singular/svn/trunk@12319 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 71.9 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009/11/03 06:57:32 $
5$Header: /usr/local/Singular/cvsroot/kernel/gfan.cc,v 1.103 2009/11/03 06:57:32 monerjan Exp $
6$Id$
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_GFAN
12
13#include "kstd1.h"
14#include "kutil.h"
15#include "intvec.h"
16#include "polys.h"
17#include "ideals.h"
18#include "kmatrix.h"
19#include "fast_maps.h"  //Mapping of ideals
20#include "maps.h"
21#include "ring.h"
22#include "structs.h"
23#include "../Singular/lists.h"
24#include "prCopy.h"
25#include <iostream>
26#include <bitset>
27#include <fstream>      //read-write cones to files
28#include <gmp.h>
29#include <string>
30#include <sstream>
31#include <time.h>
32
33//#include <gmpxx.h>
34
35/*DO NOT REMOVE THIS*/
36#ifndef GMPRATIONAL
37#define GMPRATIONAL
38#endif
39
40//Hacks for different working places
41#define p800
42
43#ifdef UNI
44#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
45#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
46#endif
47
48#ifdef HOME
49#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
50#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
51#endif
52
53#ifdef ITWM
54#include "/u/slg/monerjan/cddlib/include/setoper.h"
55#include "/u/slg/monerjan/cddlib/include/cdd.h"
56#include "/u/slg/monerjan/cddlib/include/cddmp.h"
57#endif
58
59#ifdef p800
60#include "../../cddlib/include/setoper.h"
61#include "../../cddlib/include/cdd.h"
62#include "../../cddlib/include/cddmp.h"
63#endif
64
65#ifndef gfan_DEBUG
66// #define gfan_DEBUG 1
67#ifndef gfan_DEBUGLEVEL
68#define gfan_DEBUGLEVEL 1
69#endif
70#endif
71
72#include <gfan.h>
73using namespace std;
74
75/**
76*\brief Class facet
77*       Implements the facet structure as a linked list
78*
79*/
80               
81/** \brief The default constructor for facets
82*/
83facet::facet()                 
84{
85        // Pointer to next facet.  */
86        /* Defaults to NULL. This way there is no need to check explicitly */
87        this->fNormal=NULL;
88        this->interiorPoint=NULL;               
89        this->UCN=0;
90        this->codim2Ptr=NULL;
91        this->codim=1;          //default for (codim-1)-facets
92        this->numCodim2Facets=0;
93        this->flipGB=NULL;
94        this->isIncoming=FALSE;
95        this->next=NULL;
96        this->prev=NULL;               
97        this->flipRing=NULL;    //the ring on the other side
98        this->isFlippable=FALSE;
99}
100               
101/** \brief Constructor for facets of codim >= 2
102* Note that as of now the code of the constructors is only for facets and codim2-faces. One
103* could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
104*/
105facet::facet(int const &n)
106{
107        this->fNormal=NULL;
108        this->interiorPoint=NULL;                       
109        this->UCN=0;
110        this->codim2Ptr=NULL;
111        if(n>1)
112        {
113                this->codim=n;
114        }//NOTE Handle exception here!                 
115        this->numCodim2Facets=0;
116        this->flipGB=NULL;
117        this->isIncoming=FALSE;
118        this->next=NULL;
119        this->prev=NULL;
120        this->flipRing=NULL;
121        this->isFlippable=FALSE;
122}
123               
124/** \brief The copy constructor
125*/
126facet::facet(const facet& f)
127{
128//      this->setFacetNormal(f.getFacetNormal());
129//      this->setUCN(f.getUCN());
130//      this->isFlippable=f.isFlippable;
131//      facet* f2Copy;
132//      f2Copy=f.codim2Ptr;
133//      facet* f2Act;
134//      f2Act=this->codim2Ptr;
135//      while(f2Copy!=NULL)
136//      {
137//              if(f2Act==NULL)
138//              {
139//                      f2Act=new facet(2);                     
140//              }
141//              else
142//              {
143//                      f2Act->next = new facet(2);
144//              }
145//              f2Act->setFacetNormal(f2Copy->getFacetNormal());
146//              f2Copy = f2Copy->next;
147//      }
148}
149               
150/** The default destructor */
151facet::~facet()
152{
153        delete this->fNormal;
154        delete this->interiorPoint;
155        /* Cleanup the codim2-structure */
156        if(this->codim==2)
157        {
158                facet *codim2Ptr;
159                codim2Ptr = this->codim2Ptr;
160                while(codim2Ptr!=NULL)
161                {
162                        delete codim2Ptr->fNormal;
163                        codim2Ptr = codim2Ptr->next;
164                }
165        }
166        delete this->codim2Ptr;
167        idDelete((ideal *)&this->flipGB);
168        rDelete(this->flipRing);
169        this->flipRing=NULL;
170        //this=NULL;
171}
172               
173               
174/** \brief Comparison of facets*/
175bool facet::areEqual(facet &f, facet &g)
176{
177        bool res = TRUE;
178        intvec *fNormal = new intvec(pVariables);
179        intvec *gNormal = new intvec(pVariables);
180        fNormal = f.getFacetNormal();
181        gNormal = g.getFacetNormal();
182        if((fNormal == gNormal))//||(gcone::isParallel(fNormal,gNormal)))
183        {
184                if(f.numCodim2Facets==g.numCodim2Facets)
185                {
186                        facet *f2Act;
187                        facet *g2Act;
188                        f2Act = f.codim2Ptr;
189                        g2Act = g.codim2Ptr;
190                        intvec *f2N = new intvec(pVariables);
191                        intvec *g2N = new intvec(pVariables);
192                        while(f2Act->next!=NULL && g2Act->next!=NULL)
193                        {
194                                for(int ii=0;ii<f2N->length();ii++)
195                                {
196                                        if(f2Act->getFacetNormal() != g2Act->getFacetNormal())
197                                        {
198                                                res=FALSE;                                                             
199                                        }
200                                        if (res==FALSE)
201                                                break;
202                                }
203                                if(res==FALSE)
204                                        break;
205                               
206                                f2Act = f2Act->next;
207                                g2Act = g2Act->next;
208                        }//while
209                }//if           
210        }//if
211        else
212        {
213                res = FALSE;
214        }
215        return res;
216}
217               
218/** Stores the facet normal \param intvec*/
219void facet::setFacetNormal(intvec *iv)
220{
221        this->fNormal = ivCopy(iv);                     
222}
223               
224/** Hopefully returns the facet normal */
225intvec *facet::getFacetNormal()
226{                               
227        return ivCopy(this->fNormal);
228}
229               
230/** Method to print the facet normal*/
231void facet::printNormal()
232{
233        fNormal->show();
234}
235               
236/** Store the flipped GB*/
237void facet::setFlipGB(ideal I)
238{
239        this->flipGB=idCopy(I);
240}
241               
242/** Return the flipped GB
243Seems not be used
244Anyhow idCopy would make sense here.
245*/
246ideal facet::getFlipGB()
247{
248        return this->flipGB;
249}
250               
251/** Print the flipped GB*/
252void facet::printFlipGB()
253{
254#ifndef NDEBUG
255        idShow(this->flipGB);
256#endif
257}
258               
259/** Set the UCN */
260void facet::setUCN(int n)
261{
262        this->UCN=n;
263}
264               
265/** \brief Get the UCN
266* Returns the UCN iff this != NULL, else -1
267*/
268int facet::getUCN()
269{
270        if(this!=NULL  && this!=(facet *)0xfbfbfbfbfbfbfbfb)// || this!=(facet *)0xfbfbfbfb) )
271//      if(this!=NULL && ( this->fNormal!=(intvec *)0xfbfbfbfb || this->fNormal!=(intvec *)0xfbfbfbfbfbfbfbfb) )
272                return this->UCN;
273        else
274                return -1;
275}
276               
277/** Store an interior point of the facet */
278void facet::setInteriorPoint(intvec *iv)
279{
280        this->interiorPoint = ivCopy(iv);
281}
282               
283intvec *facet::getInteriorPoint()
284{
285        return this->interiorPoint;
286}       
287               
288/** \brief Debugging function
289* prints the facet normal an all (codim-2)-facets that belong to it
290*/
291void facet::fDebugPrint()
292{                       
293        facet *codim2Act;                       
294        codim2Act = this->codim2Ptr;
295        intvec *fNormal;
296        fNormal = this->getFacetNormal();
297        intvec *f2Normal;
298        cout << "=======================" << endl;
299        cout << "Facet normal = (";
300        fNormal->show(1,1);
301        cout << ")"<<endl;     
302        cout << "-----------------------" << endl;
303        cout << "Codim2 facets:" << endl;
304        while(codim2Act!=NULL)
305        {
306                f2Normal = codim2Act->getFacetNormal();
307                cout << "(";
308                f2Normal->show(1,0);
309                cout << ")" << endl;
310                codim2Act = codim2Act->next;
311        }
312        cout << "=======================" << endl;
313}               
314               
315//friend class gcone;   //Bad style
316
317
318/**
319*\brief Implements the cone structure
320*
321* A cone is represented by a linked list of facet normals
322* @see facet
323*/
324
325
326/** \brief Default constructor.
327 *
328 * Initialises this->next=NULL and this->facetPtr=NULL
329 */
330gcone::gcone()
331{
332        this->next=NULL;
333        this->prev=NULL;
334        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();                   
335        this->baseRing=currRing;
336        this->counter++;
337        this->UCN=this->counter;                       
338        this->numFacets=0;
339        this->ivIntPt=NULL;
340}
341               
342/** \brief Constructor with ring and ideal
343 *
344 * This constructor takes the root ring and the root ideal as parameters and stores
345 * them in the private members gcone::rootRing and gcone::inputIdeal
346 * This constructor is only called once in the computation of the Gröbner fan,
347 * namely for the very first cone. Therefore pred is set to 1.
348 * Might set it to this->UCN though...
349 * Since knowledge of the root ring is only needed when using reverse search,
350 * this constructor is not needed when using the "second" method
351*/
352gcone::gcone(ring r, ideal I)
353{
354        this->next=NULL;
355        this->prev=NULL;
356        this->facetPtr=NULL;                   
357//      this->rootRing=r;
358        this->inputIdeal=I;
359        this->baseRing=currRing;
360        this->counter++;
361        this->UCN=this->counter;
362        this->pred=1;
363        this->numFacets=0;
364        this->ivIntPt=NULL;
365}
366               
367/** \brief Copy constructor
368 *
369 * Copies a cone, sets this->gcBasis to the flipped GB
370 * Call this only after a successful call to gcone::flip which sets facet::flipGB
371*/             
372gcone::gcone(const gcone& gc, const facet &f)
373{
374        this->next=NULL;
375//      this->prev=(gcone *)&gc; //comment in to get a tree
376        this->prev=NULL;
377        this->numVars=gc.numVars;                                               
378        this->counter++;
379        this->UCN=this->counter;
380        this->pred=gc.UCN;                     
381        this->facetPtr=NULL;
382        this->gcBasis=idCopy(f.flipGB);
383        this->inputIdeal=idCopy(this->gcBasis);
384        this->baseRing=rCopy(f.flipRing);
385        this->numFacets=0;
386        this->ivIntPt=NULL;
387//      this->rootRing=NULL;
388        //rComplete(this->baseRing);
389        //rChangeCurrRing(this->baseRing);
390}
391               
392/** \brief Default destructor
393*/
394gcone::~gcone()
395{
396        if(this->gcBasis!=NULL)
397                idDelete((ideal *)&this->gcBasis);
398        if(this->inputIdeal!=NULL)
399                idDelete((ideal *)&this->inputIdeal);
400//      if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
401//              rDelete(this->rootRing);
402//      if(this->UCN!=1)
403//              rDelete(this->baseRing);
404        facet *fAct;
405        facet *fDel;
406        /*Delete the facet structure*/
407        fAct=this->facetPtr;
408        fDel=fAct;
409        while(fAct!=NULL)
410        {
411                fDel=fAct;
412                fAct=fAct->next;
413                delete fDel;
414        }
415        this->counter--;
416        dd_FreeMatrix(this->ddFacets);
417        //dd_FreeMatrix(this->ddFacets);
418}                       
419       
420int gcone::getCounter()
421{
422        return this->counter;
423}
424       
425/** \brief Set the interior point of a cone */
426void gcone::setIntPoint(intvec *iv)
427{
428        this->ivIntPt=ivCopy(iv);
429}
430               
431/** \brief Return the interior point */
432intvec *gcone::getIntPoint()
433{
434        return ivCopy(this->ivIntPt);
435}
436               
437/** \brief Print the interior point */
438void gcone::showIntPoint()
439{
440        ivIntPt->show();
441}
442               
443/** \brief Print facets
444 * This is mainly for debugging purposes. Usually called from within gdb
445 */
446void gcone::showFacets(short codim)
447{
448        facet *f;
449        switch(codim)
450        {
451                case 1:
452                        f = this->facetPtr;
453                        break;
454                case 2:
455                        f = this->facetPtr->codim2Ptr;
456                        break;
457        }                       
458                       
459        intvec *iv;                     
460        while(f!=NULL)
461        {
462                iv = f->getFacetNormal();
463                cout << "(";
464                iv->show(1,0);                         
465                if(f->isFlippable==FALSE)
466                        cout << ")* ";
467                else
468                        cout << ") ";
469                f=f->next;
470        }
471        cout << endl;
472}
473               
474/** For debugging purposes only */
475volatile void gcone::showSLA(facet &f)
476{
477        facet *fAct;
478        fAct = &f;
479        facet *codim2Act;
480        codim2Act = fAct->codim2Ptr;
481        intvec *fNormal;
482        intvec *f2Normal;
483        cout << endl;
484        while(fAct!=NULL)
485        {
486                fNormal=fAct->getFacetNormal();
487                cout << "(";
488                fNormal->show(1,0);
489                if(fAct->isFlippable==TRUE)
490                        cout << ") ";
491                else
492                        cout << ")* ";
493                codim2Act = fAct->codim2Ptr;
494                cout << " Codim2: ";
495                while(codim2Act!=NULL)
496                {
497                        f2Normal = codim2Act->getFacetNormal();
498                        cout << "(";
499                        f2Normal->show(1,0);
500                        cout << ") ";
501                        codim2Act = codim2Act->next;
502                }
503                cout << "UCN = " << fAct->getUCN() << endl;                             
504                fAct = fAct->next;
505        }
506}
507               
508void gcone::idDebugPrint(ideal const &I)
509{
510        int numElts=IDELEMS(I);
511        cout << "Ideal with " << numElts << " generators" << endl;
512        cout << "Leading terms: ";
513        for (int ii=0;ii<numElts;ii++)
514        {
515                pWrite0(pHead(I->m[ii]));
516                cout << ",";
517        }
518        cout << endl;
519}
520
521void gcone::invPrint(ideal const &I)
522{
523//      int numElts=IDELEMS(I);
524//      cout << "inv = ";
525//      for(int ii=0;ii<numElts;ii++);
526//      {
527//              pWrite0(pHead(I->m[ii]));
528//              cout << ",";
529//      }
530//      cout << endl;
531}
532
533bool gcone::isMonomial(ideal const &I)
534{
535        bool res = TRUE;
536        for(int ii=0;ii<IDELEMS(I);ii++)
537        {
538                if(pLength((poly)I->m[ii])>1)
539                {
540                        res = FALSE;
541                        break;
542                }                                               
543        }
544        return res;
545}
546               
547/** \brief Set gcone::numFacets */
548void gcone::setNumFacets()
549{
550}
551               
552/** \brief Get gcone::numFacets */
553int gcone::getNumFacets()
554{
555        return this->numFacets;
556}
557               
558int gcone::getUCN()
559{
560        if( this!=NULL && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
561                return this->UCN;
562        else
563                return -1;
564}
565
566int gcone::getPredUCN()
567{
568        return this->pred;
569}
570
571/** \brief Compute the normals of the cone
572 *
573 * This method computes a representation of the cone in terms of facet normals. It takes an ideal
574 * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
575 * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
576 * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
577 * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
578 * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
579 * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
580 *
581 * Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
582 * an interior point of the cone.
583                 */
584void gcone::getConeNormals(ideal const &I, bool compIntPoint)
585{
586        poly aktpoly;
587        int rows;                       // will contain the dimensions of the ineq matrix - deprecated by
588        dd_rowrange ddrows;
589        dd_colrange ddcols;
590        dd_rowset ddredrows;            // # of redundant rows in ddineq
591        dd_rowset ddlinset;             // the opposite
592        dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
593        dd_NumberType ddnumb=dd_Integer; //Number type
594        dd_ErrorType dderr=dd_NoError;         
595        //Compute the # inequalities i.e. rows of the matrix
596        rows=0; //Initialization
597        for (int ii=0;ii<IDELEMS(I);ii++)
598        {
599                aktpoly=(poly)I->m[ii];
600                rows=rows+pLength(aktpoly)-1;
601        }
602
603        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
604        dd_set_global_constants();
605        ddrows=rows;
606        ddcols=this->numVars;
607        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
608        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
609               
610        // We loop through each g\in GB and compute the resulting inequalities
611        for (int i=0; i<IDELEMS(I); i++)
612        {
613                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
614                //simpler version of storing expvect diffs
615                int *leadexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
616                int *tailexpv=(int*)omAlloc(((this->numVars)+1)*sizeof(int));
617                pGetExpV(aktpoly,leadexpv);
618                while(pNext(aktpoly)!=NULL)
619                {
620                        aktpoly=pNext(aktpoly);
621                        pGetExpV(aktpoly,tailexpv);
622                        for(int kk=1;kk<=this->numVars;kk++)
623                        {
624                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
625                        }
626                        aktmatrixrow += 1;
627                }
628                omFree(tailexpv);
629                omFree(leadexpv);       
630        } //for
631#if true
632        /*Let's make the preprocessing here. This could already be done in the above for-loop,
633        * but for a start it is more convenient here.
634        * We check the necessary condition of FJT p.18
635        * Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
636        */
637        ideal initialForm=idInit(IDELEMS(I),1);
638        intvec *gamma=new intvec(this->numVars);
639        int falseGammaCounter=0;
640        int *redRowsArray=NULL;
641        int num_alloc=0;
642        int num_elts=0; 
643        for(int ii=0;ii<ddineq->rowsize;ii++)
644        {
645                //read row ii into gamma
646                double tmp;     
647                for(int jj=1;jj<=this->numVars;jj++)
648                {
649                        tmp=mpq_get_d(ddineq->matrix[ii][jj]);
650                        (*gamma)[jj-1]=(int)tmp;
651                }
652                computeInv((ideal&)I,initialForm,*gamma);
653                //Create leading ideal
654                ideal L=idInit(IDELEMS(initialForm),1);
655                for(int jj=0;jj<IDELEMS(initialForm);jj++)
656                {
657                        L->m[jj]=pCopy(pHead(initialForm->m[jj]));
658                }               
659               
660                LObject *P = new sLObject();
661                memset(P,0,sizeof(LObject));
662
663                for(int jj=0;jj<=IDELEMS(initialForm)-2;jj++)
664                {
665                        bool isMaybeFacet=FALSE;
666                        P->p1=initialForm->m[jj];       //build spolys of initialForm in_v
667
668                        for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
669                        {
670                                P->p2=initialForm->m[kk];
671                                ksCreateSpoly(P);                                                       
672                                if(P->p!=NULL)  //spoly non zero=?
673                                {       
674                                        poly p=pInit();                 
675                                        poly q=pInit();
676                                        p=pCopy(P->p);
677                                        q=pHead(p);     //Monomial q
678                                        isMaybeFacet=FALSE;
679                                        while(p!=NULL)
680                                        {
681                                                q=pHead(p);                                             
682//                                              unsigned long sevSpoly=pGetShortExpVector(q);
683//                                              unsigned long not_sevL;                                         
684                                                for(int ll=0;ll<IDELEMS(L);ll++)
685                                                {
686//                                                      not_sevL=~pGetShortExpVector(L->m[ll]);//                                       
687                                                        //if(!(sevSpoly & not_sevL) && pLmDivisibleBy(L->m[ll],q) )//i.e. spoly is in L
688                                                        if(pLmEqual(L->m[ll],q) || pDivisibleBy(L->m[ll],q))
689                                                        {                                                       
690                                                                isMaybeFacet=TRUE;
691                                                                break;//for
692                                                        }
693                                                }
694                                                if(isMaybeFacet==TRUE)
695                                                {                                                       
696                                                        break;//while(p!=NULL)
697                                                }
698                                                p=pNext(p);                                             
699                                        }//while
700                                        if(isMaybeFacet==FALSE)
701                                        {
702                                                dd_set_si(ddineq->matrix[ii][0],1);                                             
703                                                if(num_alloc==0)
704                                                        num_alloc += 1;
705                                                else
706                                                        num_alloc += 1;
707                                                void *_tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
708                                                if(!_tmp)
709                                                        WerrorS("Woah dude! Couldn't realloc memory\n");
710                                                redRowsArray = (int*)_tmp;
711                                                redRowsArray[num_elts]=ii;
712                                                num_elts++;
713                                                //break;//for(int kk, since we have found one that is not in L 
714                                                goto _start;    //mea culpa, mea culpa, mea maxima culpa
715                                        }
716                                }//if(P->p!=NULL)                               
717                        }//for k
718                }//for jj
719                _start:;
720                delete P;
721                //idDelete(L);         
722        }//for(ii<ddineq-rowsize
723       
724        int offset=0;//needed for correction of redRowsArray[ii]
725        for( int ii=0;ii<num_elts;ii++ )
726        {
727                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
728                offset++;
729        }
730        free(redRowsArray);
731#endif
732//      cout << "Found " << falseGammaCounter << " false in " << rows << endl;
733        //Maybe add another row to contain the constraints of the standard simplex?
734
735
736        //Remove redundant rows here!
737        //Necessary check here! C.f. FJT p18           
738        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
739        //time(&canonicalizeTac);
740        //cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl;
741        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
742        ddcols = ddineq->colsize;
743                       
744        //ddCreateMatrix(ddrows,ddcols+1);
745        ddFacets = dd_CopyMatrix(ddineq);
746                       
747        /*Write the normals into class facet*/
748        /*The pointer *fRoot should be the return value of this function*/
749        //facet *fRoot = new facet();   //instantiate new facet
750        //fRoot->setUCN(this->getUCN());
751        //this->facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
752        facet *fAct;    //pointer to active facet
753                        //fAct = fRoot;                 //Seems to do the trick. fRoot and fAct have to point to the same adress!
754                        //std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
755        int numNonFlip=0;
756        for (int kk = 0; kk<ddrows; kk++)
757        {
758                intvec *load = new intvec(this->numVars);       //intvec to store a single facet normal that will then be stored via setFacetNormal
759                for (int jj = 1; jj <ddcols; jj++)
760                {
761                        double foo;
762                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
763                        (*load)[jj-1] = (int)foo;       //store typecasted entry at pos jj-1 of load           
764                }//for (int jj = 1; jj <ddcols; jj++)
765                               
766                /*Quick'n'dirty hack for flippability*/ 
767                bool isFlip=FALSE;                                                             
768                for (int jj = 0; jj<load->length(); jj++)
769                {
770                        intvec *ivCanonical = new intvec(load->length());
771                        (*ivCanonical)[jj]=1;
772//                                      cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl;
773                        if (dotProduct(*load,*ivCanonical)<0)   
774                                        //if (ivMult(load,ivCanonical)<0)
775                        {
776                                isFlip=TRUE;
777                                break;  //URGHS
778                        }
779                        delete ivCanonical;
780                }/*End of check for flippability*/
781                if (isFlip==FALSE)
782                {
783                        this->numFacets++;
784                        numNonFlip++;
785                        if(this->numFacets==1)
786                        {
787                                facet *fRoot = new facet();
788                                this->facetPtr = fRoot;
789                                fAct = fRoot;
790                                               
791                        }
792                        else
793                        {
794                                fAct->next = new facet();
795                                fAct = fAct->next;
796                        }
797                        fAct->isFlippable=FALSE;
798                        fAct->setFacetNormal(load);
799                        fAct->setUCN(this->getUCN());
800#ifdef gfan_DEBUG
801                        cout << "Marking facet (";
802                        load->show(1,0);
803                        cout << ") as non flippable" << endl;                           
804#endif
805                }
806                else
807                {       /*Now load should be full and we can call setFacetNormal*/                                     
808                        this->numFacets++;
809                        if(this->numFacets==1)
810                        {
811                                facet *fRoot = new facet();
812                                this->facetPtr = fRoot;
813                                fAct = fRoot;
814                        }
815                        else
816                        {
817                                fAct->next = new facet();
818                                fAct = fAct->next;
819                        }
820                        fAct->isFlippable=TRUE;
821                        fAct->setFacetNormal(load);
822                        fAct->setUCN(this->getUCN());                                   
823                }//if (isFlippable==FALSE)
824                delete load;                           
825        }//for (int kk = 0; kk<ddrows; kk++)
826                       
827        //In cases like I=<x-1,y-1> there are only non-flippable facets...
828        if(numNonFlip==this->numFacets)
829        {                                       
830                WerrorS ("Only non-flippable facets. Terminating...\n");
831        }
832                       
833        /*
834        Now we should have a linked list containing the facet normals of those facets that are
835        -irredundant
836        -flipable
837        Adressing is done via *facetPtr
838        */
839                       
840        if (compIntPoint==TRUE)
841        {
842                intvec *iv = new intvec(this->numVars);
843                dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
844                int jj=1;
845                for (int ii=0;ii<=this->numVars;ii++)
846                {
847                        dd_set_si(posRestr->matrix[ii][jj],1);
848                        jj++;                                                   
849                }
850                dd_MatrixAppendTo(&ddineq,posRestr);
851                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
852                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
853                delete iv;
854        }
855                       
856        //Compute the number of facets
857        // wrong because ddineq->rowsize contains only irredundant facets. There can still be
858        // non-flippable ones!
859        //this->numFacets=ddineq->rowsize;
860       
861        //Clean up but don't delete the return value! (Whatever it will turn out to be)                 
862        dd_FreeMatrix(ddineq);
863        set_free(ddredrows);
864        set_free(ddlinset);
865        //free(ddnewpos);
866        //set_free(ddlinset);
867        //NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
868        //THIS SUCKS
869        //dd_free_global_constants();
870
871}//method getConeNormals(ideal I)       
872               
873/** \brief Compute the (codim-2)-facets of a given cone
874 * This method is used during noRevS
875 * Additionally we check whether the codim2-facet normal is strictly positive. Otherwise
876 * the facet is marked as non-flippable.
877 */
878void gcone::getCodim2Normals(gcone const &gc)
879{
880        //this->facetPtr->codim2Ptr = new facet(2);     //instantiate a (codim-2)-facet
881        facet *fAct;
882        fAct = this->facetPtr;         
883        facet *codim2Act;
884        //codim2Act = this->facetPtr->codim2Ptr;
885                       
886        dd_set_global_constants();
887                       
888        dd_MatrixPtr ddineq,P,ddakt;
889        dd_rowset impl_linset, redset;
890        dd_ErrorType err;
891        dd_rowindex newpos;             
892
893        //ddineq = facets2Matrix(gc);   //get a matrix representation of the cone
894        ddineq = dd_CopyMatrix(gc.ddFacets);
895                               
896        /*Now set appropriate linearity*/
897        dd_PolyhedraPtr ddpolyh;
898        for (int ii=0; ii<this->numFacets; ii++)                       
899        {                               
900                ddakt = dd_CopyMatrix(ddineq);
901                ddakt->representation=dd_Inequality;
902                set_addelem(ddakt->linset,ii+1);                               
903                dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                   
904                ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
905                P=dd_CopyGenerators(ddpolyh);                           
906                /* We loop through each row of P normalize it by making all
907                * entries integer ones and add the resulting vector to the
908                * int matrix facet::codim2Facets */
909                for (int jj=1;jj<=P->rowsize;jj++)
910                {                                       
911                        fAct->numCodim2Facets++;
912                        if(fAct->numCodim2Facets==1)                                   
913                        {                                               
914                                fAct->codim2Ptr = new facet(2);                                         
915                                codim2Act = fAct->codim2Ptr;
916                        }
917                        else
918                        {
919                                codim2Act->next = new facet(2);
920                                codim2Act = codim2Act->next;
921                        }
922                        intvec *n = new intvec(this->numVars);
923                        makeInt(P,jj,*n);
924                        codim2Act->setFacetNormal(n);
925                        delete n;                                       
926                }               
927                /*We check whether the facet spanned by the codim-2 facets
928                * intersects with the positive orthant. Otherwise we define this
929                * facet to be non-flippable
930                */
931                intvec *iv_intPoint = new intvec(this->numVars);
932                dd_MatrixPtr shiftMatrix;
933                dd_MatrixPtr intPointMatrix;
934                shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
935                for(int kk=0;kk<this->numVars;kk++)
936                {
937                        dd_set_si(shiftMatrix->matrix[kk][0],1);
938                        dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
939                }
940                intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);             
941                interiorPoint(intPointMatrix,*iv_intPoint);
942                for(int ll=0;ll<this->numVars;ll++)
943                {
944                        if( (*iv_intPoint)[ll] < 0 )
945                        {
946                                fAct->isFlippable=FALSE;
947                                break;
948                        }
949                }
950                /*End of check*/                                       
951                fAct = fAct->next;     
952                dd_FreeMatrix(ddakt);
953//              dd_FreeMatrix(ddineq);
954                dd_FreeMatrix(shiftMatrix);
955                dd_FreeMatrix(intPointMatrix);
956                dd_FreePolyhedra(ddpolyh);
957                delete iv_intPoint;
958        }//for
959        dd_FreeMatrix(ddineq);
960}
961               
962/** \brief Compute the Groebner Basis on the other side of a shared facet
963 *
964 * Implements algorithm 4.3.2 from Anders' thesis.
965 * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
966 * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
967 * is parallel to \f$ leadexp(g) \f$
968 * Parallelity is checked using basic linear algebra. See gcone::isParallel.
969 * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
970 * computing an interior point of the facet and taking all terms having the same weight with respect
971 * to this interior point.
972 *\param ideal, facet
973 * Input: a marked,reduced Groebner basis and a facet
974 */
975void gcone::flip(ideal gb, facet *f)            //Compute "the other side"
976{                       
977        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
978        fNormal = f->getFacetNormal();  //read this->fNormal;
979
980//      std::cout << "running gcone::flip" << std::endl;
981//      std::cout << "flipping UCN " << this->getUCN() << endl;
982//      cout << "over facet (";
983//      fNormal->show(1,0);
984//      cout << ") with UCN " << f->getUCN();
985//      std::cout << std::endl;
986
987        /*1st step: Compute the initial ideal*/
988        /*poly initialFormElement[IDELEMS(gb)];*/       //array of #polys in GB to store initial form
989        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
990        poly aktpoly;   
991        computeInv(gb,initialForm,*fNormal);
992        //NOTE The code below went into gcone::computeInv
993#ifdef gfan_DEBUG
994/*      cout << "Initial ideal is: " << endl;
995        idShow(initialForm);
996        //f->printFlipGB();*/
997//      cout << "===" << endl;
998#endif                 
999        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
1000        /*Substep 2.1
1001        compute $G_{-\alpha}(in_v(I))
1002        see journal p. 66
1003        NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
1004        srcRing already has a weighted ordering
1005        */
1006        ring srcRing=currRing;
1007        ring tmpRing;
1008                       
1009        if( (srcRing->order[0]!=ringorder_a))
1010        {
1011                tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
1012        }
1013        else
1014        {
1015                tmpRing=rCopy0(srcRing);
1016                int length=fNormal->length();
1017                int *A=(int *)omAlloc0(length*sizeof(int));
1018                for(int jj=0;jj<length;jj++)
1019                {
1020                        A[jj]=-(*fNormal)[jj];
1021                }
1022                omFree(tmpRing->wvhdl[0]);
1023                tmpRing->wvhdl[0]=(int*)A;
1024                tmpRing->block1[0]=length;
1025                rComplete(tmpRing);     
1026        }
1027        delete fNormal; //NOTE: Why does this crash?
1028        rChangeCurrRing(tmpRing);
1029                       
1030        //rWrite(currRing); cout << endl;
1031                       
1032        ideal ina;                     
1033        ina=idrCopyR(initialForm,srcRing);
1034#ifndef NDEBUG                 
1035#ifdef gfan_DEBUG
1036        cout << "ina=";
1037        idShow(ina); cout << endl;
1038#endif
1039#endif
1040        ideal H;
1041        //H=kStd(ina,NULL,isHomog,NULL);        //we know it is homogeneous
1042        H=kStd(ina,NULL,testHomog,NULL);        //This is \mathcal(G)_{>_-\alpha}(in_v(I))
1043        idSkipZeroes(H);
1044#ifndef NDEBUG
1045#ifdef gfan_DEBUG
1046                        cout << "H="; idShow(H); cout << endl;
1047#endif
1048#endif
1049        /*Substep 2.2
1050        do the lifting and mark according to H
1051        */
1052        rChangeCurrRing(srcRing);
1053        ideal srcRing_H;
1054        ideal srcRing_HH;                       
1055        srcRing_H=idrCopyR(H,tmpRing);
1056#ifndef NDEBUG
1057#ifdef gfan_DEBUG
1058        cout << "srcRing_H = ";
1059        idShow(srcRing_H); cout << endl;
1060#endif
1061#endif
1062        srcRing_HH=ffG(srcRing_H,this->gcBasis);       
1063#ifndef NDEBUG
1064#ifdef gfan_DEBUG
1065        cout << "srcRing_HH = ";
1066        idShow(srcRing_HH); cout << endl;
1067#endif
1068#endif
1069        /*Substep 2.2.1
1070         * Mark according to G_-\alpha
1071         * Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
1072         * we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
1073         * represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
1074         * Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
1075         * compute the difference accordingly
1076        */
1077        dd_set_global_constants();
1078        bool markingsAreCorrect=FALSE;
1079        dd_MatrixPtr intPointMatrix;
1080        int iPMatrixRows=0;
1081        dd_rowrange aktrow=0;                   
1082        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1083        {
1084                poly aktpoly=(poly)srcRing_HH->m[ii];
1085                iPMatrixRows = iPMatrixRows+pLength(aktpoly);
1086        }
1087        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
1088         * construction of the differences
1089         */
1090        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); //iPMatrixRows+10;
1091        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
1092       
1093        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1094        {
1095                markingsAreCorrect=FALSE;       //crucial to initialise here
1096                poly aktpoly=srcRing_HH->m[ii];
1097                /*Comparison of leading monomials is done via exponent vectors*/
1098                for (int jj=0;jj<IDELEMS(H);jj++)
1099                {
1100                        int *src_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
1101                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
1102                        pGetExpV(aktpoly,src_ExpV);
1103                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
1104                        pGetExpV(pCopy(H->m[ii]),dst_ExpV);
1105                        rChangeCurrRing(srcRing);
1106                        bool expVAreEqual=TRUE;
1107                        for (int kk=1;kk<=this->numVars;kk++)
1108                        {
1109#ifdef gfan_DEBUG
1110//                                              cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
1111#endif
1112                                if (src_ExpV[kk]!=dst_ExpV[kk])
1113                                {
1114                                        expVAreEqual=FALSE;
1115                                }
1116                        }                                       
1117                                        //if (*src_ExpV == *dst_ExpV)
1118                        if (expVAreEqual==TRUE)
1119                        {
1120                                markingsAreCorrect=TRUE; //everything is fine
1121#ifdef gfan_DEBUG
1122//                              cout << "correct markings" << endl;
1123#endif
1124                        }//if (pHead(aktpoly)==pHead(H->m[jj])
1125                        omFree(src_ExpV);
1126                        omFree(dst_ExpV);
1127                }//for (int jj=0;jj<IDELEMS(H);jj++)
1128                               
1129                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
1130                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
1131                if (markingsAreCorrect==TRUE)
1132                {
1133                        pGetExpV(aktpoly,leadExpV);
1134                }
1135                else
1136                {
1137                        rChangeCurrRing(tmpRing);
1138                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
1139                        rChangeCurrRing(srcRing);
1140                }
1141                /*compute differences of the expvects*/                         
1142                while (pNext(aktpoly)!=NULL)
1143                {
1144                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
1145                        is not omitted when computing the differences*/
1146                        if(markingsAreCorrect==TRUE)
1147                        {
1148                                aktpoly=pNext(aktpoly);
1149                                pGetExpV(aktpoly,v);
1150                        }
1151                        else
1152                        {
1153                                pGetExpV(pHead(aktpoly),v);
1154                                markingsAreCorrect=TRUE;
1155                        }
1156
1157                        for (int jj=0;jj<this->numVars;jj++)
1158                        {                               
1159                                /*Store into ddMatrix*/                                                         
1160                                dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
1161                        }
1162                        aktrow +=1;
1163                }
1164                omFree(v);
1165                omFree(leadExpV);
1166        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1167        /*Now we add the constraint for the standard simplex*/
1168        /*NOTE:Might actually work without the standard simplex
1169        * No, it won't. MM
1170        */
1171        dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
1172        for (int jj=1;jj<=this->numVars;jj++)
1173        {
1174                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
1175        }
1176        //Let's make sure we compute interior points from the positive orthant
1177        dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1178        int jj=1;
1179        for (int ii=0;ii<this->numVars;ii++)
1180        {
1181                dd_set_si(posRestr->matrix[ii][jj],1);
1182                jj++;                                                   
1183        }
1184        dd_MatrixAppendTo(&intPointMatrix,posRestr);
1185        dd_FreeMatrix(posRestr);
1186#ifdef gfan_DEBUG
1187//                      dd_WriteMatrix(stdout,intPointMatrix);
1188#endif
1189        intvec *iv_weight = new intvec(this->numVars);
1190        interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
1191        dd_FreeMatrix(intPointMatrix);
1192        dd_free_global_constants();
1193                       
1194        /*Step 3
1195        turn the minimal basis into a reduced one
1196        */
1197        //ring dstRing=rCopyAndAddWeight(tmpRing,iv_weight);   
1198                       
1199        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
1200        // Thus:
1201        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
1202        ring dstRing=rCopy0(tmpRing);
1203        int length=iv_weight->length();
1204        int *A=(int *)omAlloc0(length*sizeof(int));
1205        for(int jj=0;jj<length;jj++)
1206        {
1207                A[jj]=(*iv_weight)[jj];
1208        }
1209        dstRing->wvhdl[0]=(int*)A;
1210        rComplete(dstRing);                                     
1211        rChangeCurrRing(dstRing);
1212        //rDelete(tmpRing);
1213        delete iv_weight;
1214
1215        //rWrite(dstRing); cout << endl;
1216
1217        ideal dstRing_I;                       
1218        dstRing_I=idrCopyR(srcRing_HH,srcRing);
1219        //dstRing_I=idrCopyR(inputIdeal,srcRing);
1220        //validOpts<1>=TRUE;
1221#ifdef gfan_DEBUG
1222        //idShow(dstRing_I);
1223#endif                 
1224        BITSET save=test;
1225        test|=Sy_bit(OPT_REDSB);
1226        test|=Sy_bit(OPT_REDTAIL);
1227#ifdef gfan_DEBUG
1228        test|=Sy_bit(6);        //OPT_DEBUG
1229#endif
1230        ideal tmpI;
1231        //NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
1232        //tmpI = idrCopyR(this->inputIdeal,this->baseRing);
1233        tmpI = dstRing_I;
1234        dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
1235        idNorm(dstRing_I);                     
1236                        //kInterRed(dstRing_I);
1237        idSkipZeroes(dstRing_I);
1238        test=save;
1239        /*End of step 3 - reduction*/
1240                       
1241        f->setFlipGB(dstRing_I);//store the flipped GB
1242        f->flipRing=rCopy(dstRing);     //store the ring on the other side
1243//#ifdef gfan_DEBUG
1244//      cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
1245//      this->idDebugPrint(dstRing_I);
1246//      cout << endl;
1247//#endif                       
1248        rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
1249//      rDelete(dstRing);
1250}//void flip(ideal gb, facet *f)
1251
1252/** \brief Compute initial ideal
1253 * Compute the initial ideal in_v(G) wrt a (possible) facet normal
1254 * used in gcone::getFacetNormal in order to preprocess possible facet normals
1255 * and in gcone::flip for obvious reasons.
1256*/
1257void gcone::computeInv(ideal &gb, ideal &initialForm, intvec &fNormal)
1258{
1259        intvec *check = new intvec(this->numVars);
1260        poly initialFormElement[IDELEMS(gb)];
1261        poly aktpoly;
1262       
1263        for (int ii=0;ii<IDELEMS(gb);ii++)
1264        {
1265                aktpoly = (poly)gb->m[ii];                                                             
1266                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
1267                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
1268                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
1269                initialFormElement[ii]=pHead(aktpoly);
1270                               
1271                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
1272                {
1273                        aktpoly=pNext(aktpoly); //next term
1274                        pSetm(aktpoly);
1275                        pGetExpV(aktpoly,v);           
1276                        /* Convert (int)v into (intvec)check */                 
1277                        for (int jj=0;jj<this->numVars;jj++)
1278                        {
1279                                                //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
1280                                                //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
1281                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
1282                        }
1283#ifdef gfan_DEBUG
1284//                      cout << "check=";                       
1285//                      check->show();
1286//                      cout << endl;
1287#endif
1288                        if (isParallel(*check,fNormal)) //pass *check when
1289                        {
1290//                              cout << "Parallel vector found, adding to initialFormElement" << endl;                 
1291                                initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
1292                        }                                               
1293                }//while
1294#ifdef gfan_DEBUG
1295//              cout << "Initial Form=";                               
1296//              pWrite(initialFormElement[ii]);
1297//              cout << "---" << endl;
1298#endif
1299                /*Now initialFormElement must be added to (ideal)initialForm */
1300                initialForm->m[ii]=initialFormElement[ii];
1301                omFree(leadExpV);
1302                omFree(v);             
1303        }//for
1304        delete check;
1305}
1306
1307/** \brief Compute the remainder of a polynomial by a given ideal
1308 *
1309 * Compute \f$ f^{\mathcal{G}} \f$
1310 * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
1311 * However, since we are only interested in the remainder, there is no need to
1312 * compute the factors \f$ a_i \f$
1313 */
1314//NOTE: Should be replaced by kNF or kNF2
1315//NOTE: Done
1316//NOTE: removed with r12286
1317               
1318/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
1319*/
1320//NOTE: use kNF or kNF2 instead of restOfDivision
1321ideal gcone::ffG(ideal const &H, ideal const &G)
1322{
1323//                      cout << "Entering ffG" << endl;
1324        int size=IDELEMS(H);
1325        ideal res=idInit(size,1);
1326        poly temp1, temp2, temp3;       //polys to temporarily store values for pSub
1327        for (int ii=0;ii<size;ii++)
1328        {
1329//              res->m[ii]=restOfDiv(H->m[ii],G);
1330                res->m[ii]=kNF(G, NULL,H->m[ii],0,0);
1331                temp1=H->m[ii];
1332                temp2=res->m[ii];
1333                temp3=pSub(temp1, temp2);
1334                res->m[ii]=temp3;
1335                //res->m[ii]=pSub(temp1,temp2); //buggy
1336                //pSort(res->m[ii]);
1337                //pSetm(res->m[ii]);
1338                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);       
1339        }                       
1340        return res;
1341}
1342               
1343/** \brief Compute a Groebner Basis
1344 *
1345 * Computes the Groebner basis and stores the result in gcone::gcBasis
1346 *\param ideal
1347 *\return void
1348 */
1349void gcone::getGB(ideal const &inputIdeal)             
1350{
1351        BITSET save=test;
1352        test|=Sy_bit(OPT_REDSB);
1353        test|=Sy_bit(OPT_REDTAIL);
1354        ideal gb;
1355        gb=kStd(inputIdeal,NULL,testHomog,NULL);
1356        idNorm(gb);
1357        idSkipZeroes(gb);
1358        this->gcBasis=gb;       //write the GB into gcBasis
1359        test=save;
1360}//void getGB
1361               
1362/** \brief Compute the negative of a given intvec
1363        */             
1364intvec *gcone::ivNeg(const intvec *iv)
1365{
1366        //NOTE: Can't this be done without new?
1367        intvec *res = new intvec(iv->length());
1368        res=ivCopy(iv);
1369        *res *= (int)-1;                                               
1370        return res;
1371}
1372
1373
1374/** \brief Compute the dot product of two intvecs
1375*
1376*/
1377int gcone::dotProduct(intvec const &iva, intvec const &ivb)                             
1378{                       
1379        int res=0;
1380        for (int i=0;i<this->numVars;i++)
1381        {
1382                res = res+(iva[i]*ivb[i]);
1383        }
1384        return res;
1385}//int dotProduct
1386
1387/** \brief Check whether two intvecs are parallel
1388 *
1389 * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1390 */
1391bool gcone::isParallel(intvec const &a, intvec const &b)
1392{                       
1393        int lhs,rhs;
1394        bool res;
1395        lhs=dotProduct(a,b)*dotProduct(a,b);
1396        rhs=dotProduct(a,a)*dotProduct(b,b);
1397                        //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1398        if (lhs==rhs)
1399        {
1400                res = TRUE;
1401        }
1402        else
1403        {
1404                res = FALSE;
1405        }
1406        return res;
1407}//bool isParallel
1408               
1409/** \brief Compute an interior point of a given cone
1410 * Result will be written into intvec iv.
1411 * Any rational point is automatically converted into an integer.
1412 */
1413void gcone::interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
1414{
1415        dd_LPPtr lp,lpInt;
1416        dd_ErrorType err=dd_NoError;
1417        dd_LPSolverType solver=dd_DualSimplex;
1418        dd_LPSolutionPtr lpSol=NULL;
1419        dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
1420        dd_rowindex ddnewpos;
1421        dd_NumberType numb;     
1422                        //M->representation=dd_Inequality;
1423                        //M->objective-dd_LPMin;  //Not sure whether this is needed
1424                       
1425                        //NOTE: Make this n-dimensional!
1426                        //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1427                                                                       
1428                        /*NOTE: Leave the following line commented out!
1429        * Otherwise it will cause interior points that are not strictly positive on some examples
1430        *
1431                        */
1432                        //dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1433                        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1434                        //cout << "Tick 2" << endl;
1435                        //dd_WriteMatrix(stdout,M);
1436                       
1437        lp=dd_Matrix2LP(M, &err);
1438        if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}                     
1439        if (lp==NULL){cout << "LP is NULL" << endl;}
1440#ifdef gfan_DEBUG
1441//                      dd_WriteLP(stdout,lp);
1442#endif 
1443                                       
1444        lpInt=dd_MakeLPforInteriorFinding(lp);
1445        if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1446#ifdef gfan_DEBUG
1447//                      dd_WriteLP(stdout,lpInt);
1448#endif                 
1449
1450        dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1451        if (err!=dd_NoError)
1452        {
1453                cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1454                dd_WriteErrorMessages(stdout, err);
1455        }
1456                       
1457                        //dd_LPSolve(lpInt,solver,&err);        //This will not result in a point from the relative interior
1458        if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1459                                                                       
1460                        //lpSol=dd_CopyLPSolution(lpInt);
1461        if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}                 
1462#ifdef gfan_DEBUG
1463        cout << "Interior point: ";
1464        for (int ii=1; ii<(lpSol->d)-1;ii++)
1465        {
1466                dd_WriteNumber(stdout,lpSol->sol[ii]);
1467        }
1468        cout << endl;
1469#endif
1470                       
1471        //NOTE The following strongly resembles parts of makeInt.
1472        //Maybe merge sometimes
1473        mpz_t kgV; mpz_init(kgV);
1474        mpz_set_str(kgV,"1",10);
1475        mpz_t den; mpz_init(den);
1476        mpz_t tmp; mpz_init(tmp);
1477        mpq_get_den(tmp,lpSol->sol[1]);
1478        for (int ii=1;ii<(lpSol->d)-1;ii++)
1479        {
1480                mpq_get_den(den,lpSol->sol[ii+1]);
1481                mpz_lcm(kgV,tmp,den);
1482                mpz_set(tmp, kgV);
1483        }
1484        mpq_t qkgV;
1485        mpq_init(qkgV);
1486        mpq_set_z(qkgV,kgV);
1487        for (int ii=1;ii<(lpSol->d)-1;ii++)
1488        {
1489                mpq_t product;
1490                mpq_init(product);
1491                mpq_mul(product,qkgV,lpSol->sol[ii]);
1492                iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
1493                mpq_clear(product);
1494        }
1495#ifdef gfan_DEBUG
1496//                      iv.show();
1497//                      cout << endl;
1498#endif
1499        mpq_clear(qkgV);
1500        mpz_clear(tmp);
1501        mpz_clear(den);
1502        mpz_clear(kgV);                 
1503                       
1504        dd_FreeLPSolution(lpSol);
1505        dd_FreeLPData(lpInt);
1506        dd_FreeLPData(lp);
1507        set_free(ddlinset);
1508        set_free(ddredrows);                   
1509                       
1510}//void interiorPoint(dd_MatrixPtr const &M)
1511               
1512/** \brief Copy a ring and add a weighted ordering in first place
1513 *
1514 */
1515ring gcone::rCopyAndAddWeight(ring const &r, intvec const *ivw)                         
1516{
1517        ring res=rCopy0(r);
1518        int jj;
1519                       
1520        omFree(res->order);
1521        res->order =(int *)omAlloc0(4*sizeof(int));
1522        omFree(res->block0);
1523        res->block0=(int *)omAlloc0(4*sizeof(int));
1524        omFree(res->block1);
1525        res->block1=(int *)omAlloc0(4*sizeof(int));
1526        omfree(res->wvhdl);
1527        res->wvhdl =(int **)omAlloc0(4*sizeof(int**));
1528                       
1529        res->order[0]=ringorder_a;
1530        res->block0[0]=1;
1531        res->block1[0]=res->N;;
1532        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                   
1533        res->block0[1]=1;
1534        res->block1[1]=res->N;;
1535        res->order[2]=ringorder_C;
1536
1537        int length=ivw->length();
1538        int *A=(int *)omAlloc0(length*sizeof(int));
1539        for (jj=0;jj<length;jj++)
1540        {                               
1541                A[jj]=(*ivw)[jj];                               
1542        }                       
1543        res->wvhdl[0]=(int *)A;
1544        res->block1[0]=length;
1545                       
1546        rComplete(res);
1547        return res;
1548}//rCopyAndAdd
1549               
1550ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1551{
1552        ring res=rCopy0(currRing);
1553        rComplete(res);
1554        rSetWeightVec(res,(int64*)ivw);
1555        //rChangeCurrRing(rnew);
1556        return res;
1557}
1558               
1559/** \brief Checks whether a given facet is a search facet
1560 * Determines whether a given facet of a cone is the search facet of a neighbouring cone
1561 * This is done in the following way:
1562 * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1563 * that is first crossed during the generic walk.
1564 * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1565 * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1566*/
1567//removed with r12286
1568               
1569/** \brief Check for equality of two intvecs
1570 */
1571bool gcone::areEqual(intvec const &a, intvec const &b)
1572{
1573        bool res=TRUE;
1574        for(int ii=0;ii<this->numVars;ii++)
1575        {
1576                if(a[ii]!=b[ii])
1577                {
1578                        res=FALSE;
1579                        break;
1580                }
1581        }
1582        return res;
1583}
1584               
1585/** \brief The reverse search algorithm
1586 */
1587//removed with r12286
1588               
1589/** \brief The new method of Markwig and Jensen
1590 * Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1591 * In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1592 * Keep a list of facets as a linked list containing an intvec and an integer matrix.
1593 * Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1594 * soon as we compute this facet again. Comparison of facets is done by...
1595 */
1596void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
1597{       
1598        facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
1599        facet *SearchListAct;
1600        SearchListAct = NULL;
1601        //SearchListAct = SearchListRoot;
1602                       
1603        gcone *gcAct;
1604        gcAct = &gcRoot;
1605        gcone *gcPtr;   //Pointer to end of linked list of cones
1606        gcPtr = &gcRoot;
1607        gcone *gcNext;  //Pointer to next cone to be visited
1608        gcNext = &gcRoot;
1609        gcone *gcHead;
1610        gcHead = &gcRoot;
1611                       
1612        facet *fAct;
1613        fAct = gcAct->facetPtr;                 
1614                       
1615        ring rAct;
1616        rAct = currRing;
1617                                               
1618        int UCNcounter=gcAct->getUCN();
1619       
1620        /* Use pure SLA or writeCone2File */
1621//      enum heuristic {
1622//              ram,
1623//              hdd
1624//      };
1625//      heuristic h;
1626//      h=hdd;
1627       
1628#ifdef gfan_DEBUG
1629        cout << "NoRevs" << endl;
1630        cout << "Facets are:" << endl;
1631        gcAct->showFacets();
1632#endif                 
1633                       
1634        gcAct->getCodim2Normals(*gcAct);
1635                               
1636        //Compute unique representation of codim-2-facets
1637        gcAct->normalize();
1638                       
1639        /* Make a copy of the facet list of first cone
1640           Since the operations getCodim2Normals and normalize affect the facets
1641           we must not memcpy them before these ops!
1642        */
1643        intvec *fNormal;// = new intvec(this->numVars);
1644        intvec *f2Normal;// = new intvec(this->numVars);
1645        facet *codim2Act; codim2Act = NULL;                     
1646        facet *sl2Root; //sl2Root = new facet(2);
1647        facet *sl2Act;  //sl2Act = sl2Root;
1648                       
1649        for(int ii=0;ii<this->numFacets;ii++)
1650        {
1651                //only copy flippable facets! NOTE: We do not need flipRing or any such information.
1652                if(fAct->isFlippable==TRUE)
1653                {
1654                        fNormal = fAct->getFacetNormal();
1655                        if( ii==0 || (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
1656                        {                                               
1657                                SearchListAct = SearchListRoot;
1658                                                //memcpy(SearchListAct,fAct,sizeof(facet));                                     
1659                        }
1660                        else
1661                        {                       
1662                                SearchListAct->next = new facet();
1663                                SearchListAct = SearchListAct->next;                                           
1664                                                //memcpy(SearchListAct,fAct,sizeof(facet));                             
1665                        }
1666                        SearchListAct->setFacetNormal(fNormal);
1667                        SearchListAct->setUCN(this->getUCN());
1668                        SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
1669                        SearchListAct->isFlippable=TRUE;
1670                                        //Copy codim2-facets                           
1671                        codim2Act=fAct->codim2Ptr;
1672                        SearchListAct->codim2Ptr = new facet(2);
1673                        sl2Root = SearchListAct->codim2Ptr;
1674                        sl2Act = sl2Root;
1675                                        //while(codim2Act!=NULL)
1676                        for(int jj=0;jj<fAct->numCodim2Facets;jj++)
1677                        {
1678                                f2Normal = codim2Act->getFacetNormal();
1679                                if(jj==0)
1680                                {                                               
1681                                        sl2Act = sl2Root;
1682                                        sl2Act->setFacetNormal(f2Normal);
1683                                }
1684                                else
1685                                {
1686                                        sl2Act->next = new facet(2);
1687                                        sl2Act = sl2Act->next;
1688                                        sl2Act->setFacetNormal(f2Normal);
1689                                }                                       
1690                                codim2Act = codim2Act->next;
1691                        }
1692                        fAct = fAct->next;
1693                }//if(fAct->isFlippable==TRUE)
1694                else {fAct = fAct->next;}
1695        }//End of copying facets into SLA
1696       
1697        SearchListAct = SearchListRoot; //Set to beginning of list
1698        /*Make SearchList doubly linked*/
1699        while(SearchListAct!=NULL)
1700        {
1701                if(SearchListAct->next!=NULL)
1702                {
1703                        SearchListAct->next->prev = SearchListAct;                                     
1704                }
1705                SearchListAct = SearchListAct->next;
1706        }
1707        SearchListAct = SearchListRoot; //Set to beginning of List
1708       
1709        fAct = gcAct->facetPtr;
1710        //NOTE Disabled until code works as expected. Reenabled 2.11.2009
1711        gcAct->writeConeToFile(*gcAct);
1712                       
1713        /*End of initialisation*/
1714        fAct = SearchListAct;
1715        /*2nd step
1716          Choose a facet from fListPtr, flip it and forget the previous cone
1717          We always choose the first facet from fListPtr as facet to be flipped
1718        */                     
1719//      time_t tic, tac;
1720        while((SearchListAct!=NULL))// && counter<10)
1721        {//NOTE See to it that the cone is only changed after ALL facets have been flipped!                             
1722                fAct = SearchListAct;
1723                               
1724                while(fAct!=NULL)
1725                {       //Since SLA should only contain flippables there should be no need to check for that
1726                       
1727                        gcAct->flip(gcAct->gcBasis,fAct);
1728                       
1729                       
1730                        ring rTmp=rCopy(fAct->flipRing);
1731                        rComplete(rTmp);                       
1732                        rChangeCurrRing(rTmp);
1733                        gcone *gcTmp = new gcone::gcone(*gcAct,*fAct);//copy constructor!
1734                        /* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
1735                         * Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
1736                         * NOTE: Can this cause trouble with the destructor?
1737                         * Answer: Yes it bloody well does!
1738                         * However I'd like to delete this data here, since if we have a cone with many many facets it
1739                         * should pay to get rid of the flipGB as soon as possible.
1740                         * Destructor must be equipped with necessary checks.
1741                        */
1742//                      idDelete((ideal *)&fAct->flipGB);
1743//                      rDelete(fAct->flipRing);
1744                       
1745                        gcTmp->getConeNormals(gcTmp->gcBasis, FALSE);   
1746                        gcTmp->getCodim2Normals(*gcTmp);
1747                        gcTmp->normalize();
1748#ifdef gfan_DEBUG
1749//                      gcTmp->showFacets(1);
1750#endif
1751                        /*add facets to SLA here*/
1752                        SearchListRoot=gcTmp->enqueueNewFacets(*SearchListRoot);
1753                       
1754                        if(gfanHeuristic==1)
1755                        {
1756                                gcTmp->writeConeToFile(*gcTmp);
1757//                              idDelete((ideal*)&gcTmp->gcBasis);
1758                                for(int ii=0;ii<IDELEMS(gcTmp->gcBasis);ii++)
1759                                {
1760//                                      gcTmp->gcBasis->m[ii]=(poly)NULL;
1761                                        pDelete(&gcTmp->gcBasis->m[ii]);
1762                                }
1763                        }
1764                       
1765#if gfan_DEBUG
1766//                      if(SearchListRoot!=NULL)
1767//                              gcTmp->showSLA(*SearchListRoot);
1768#endif
1769                        rChangeCurrRing(gcAct->baseRing);
1770                        //rDelete(rTmp);
1771                        //doubly linked for easier removal
1772                        gcTmp->prev = gcPtr;
1773                        gcPtr->next=gcTmp;
1774                        gcPtr=gcPtr->next;
1775                        if(fAct->getUCN() == fAct->next->getUCN())
1776                        {
1777                                fAct=fAct->next;
1778                        }
1779                        else
1780                                break;
1781                }//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
1782                //NOTE Now all neighbouring cones of gcAct have been computed, so we may delete gcAct
1783                gcone *deleteMarker;
1784                deleteMarker=gcAct;             
1785//              delete deleteMarker;
1786//              deleteMarker=NULL;
1787                //Search for cone with smallest UCN
1788                gcNext = gcHead;
1789                while(gcNext!=NULL && SearchListRoot!=NULL && gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && gcNext!=(gcone * const)0xfbfbfbfb)
1790                {                       
1791                        /*NOTE if gcNext->getUCN is smaller than SearchListRoot->getUCN it follows, that the cone
1792                        gcNext will not be used in any more computation. So -> delete
1793                        BAD IDEA!!! Luckily the delete statement below is never executed since noRevs is
1794                        invoked from gcRoot which has UCN==1. And apparently the else below is also never executed!
1795                        */
1796                        if (gcNext->getUCN() < SearchListRoot->getUCN() )
1797                        {
1798                                //idDelete((ideal *)&gcNext->gcBasis);  //at least drop huge gröbner bases
1799                                if(gcNext == gcHead)
1800                                {
1801                                        deleteMarker = gcHead;
1802//                                      gcHead = gcNext->next;
1803                                        //gcNext->next->prev = NULL;
1804                                }
1805                                else
1806                                {
1807                                        deleteMarker = gcNext;
1808//                                      gcNext->prev->next = gcNext->next;
1809//                                      gcNext->next->prev = gcNext->prev;
1810                                }
1811//                              gcNext = gcNext->next;
1812//                              delete deleteMarker;
1813//                              deleteMarker = NULL;
1814                        }
1815                        else if( gcNext->getUCN() == SearchListRoot->getUCN() )
1816                        {//NOTE: Implement heuristic to be used!
1817                                if(gfanHeuristic==0)
1818                                {
1819                                        gcAct = gcNext;
1820                                        rAct=rCopy(gcAct->baseRing);
1821                                        rComplete(rAct);
1822                                        rChangeCurrRing(rAct);                                         
1823                                        break;
1824                                }
1825                                else if(gfanHeuristic==1)
1826                                {
1827                                        //Read st00f from file
1828                                        gcAct = gcNext;
1829                                        //implant the GB into gcAct
1830                                        readConeFromFile(gcAct->getUCN(), gcAct);
1831                                        rAct=rCopy(gcAct->baseRing);
1832                                        rComplete(rAct);
1833                                        rChangeCurrRing(rAct);
1834                                        break;
1835                                }                               
1836                        }                       
1837                        gcNext = gcNext->next;
1838//                      if(deleteMarker!=NULL && deleteMarker!=(gcone *)0xfbfbfbfbfbfbfbfb)
1839                        if(this->getUCN()!=1)   //workaround to avoid unpredictable behaviour towards end of program
1840                                delete deleteMarker;
1841//                      deleteMarker = NULL;
1842//                      gcNext = gcNext->next;
1843                }
1844                UCNcounter++;
1845                                //SearchListAct = SearchListAct->next;
1846                                //SearchListAct = fAct->next;
1847                SearchListAct = SearchListRoot;
1848        }
1849        cout << endl << "Found " << counter << " cones - terminating" << endl;
1850               
1851                        //NOTE Hm, comment in and get a crash for free...
1852                        //dd_free_global_constants();                           
1853                        //gc.writeConeToFile(gc);
1854                       
1855                       
1856                        //fAct = fListPtr;
1857                        //gcone *gcTmp = new gcone(gc); //copy
1858                        //gcTmp->flip(gcTmp->gcBasis,fAct);
1859                        //NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
1860                       
1861}//void noRevS(gcone &gc)       
1862               
1863               
1864/** \brief Make a set of rational vectors into integers
1865 *
1866 * We compute the lcm of the denominators and multiply with this to get integer values.         
1867 * \param dd_MatrixPtr,intvec
1868 */
1869void gcone::makeInt(dd_MatrixPtr const &M, int const line, intvec &n)
1870{                       
1871        mpz_t denom[this->numVars];
1872        for(int ii=0;ii<this->numVars;ii++)
1873        {
1874                mpz_init_set_str(denom[ii],"0",10);
1875        }
1876               
1877        mpz_t kgV,tmp;
1878        mpz_init(kgV);
1879        mpz_init(tmp);
1880                       
1881        for (int ii=0;ii<(M->colsize)-1;ii++)
1882        {
1883                mpz_t z;
1884                mpz_init(z);
1885                mpq_get_den(z,M->matrix[line-1][ii+1]);
1886                mpz_set( denom[ii], z);
1887                mpz_clear(z);                           
1888        }
1889                                               
1890        /*Compute lcm of the denominators*/
1891        mpz_set(tmp,denom[0]);
1892        for (int ii=0;ii<(M->colsize)-1;ii++)
1893        {
1894                mpz_lcm(kgV,tmp,denom[ii]);
1895                mpz_set(tmp,kgV);                               
1896        }
1897                       
1898        /*Multiply the nominators by kgV*/
1899        mpq_t qkgV,res;
1900        mpq_init(qkgV);
1901        mpq_set_str(qkgV,"1",10);
1902        mpq_canonicalize(qkgV);
1903                       
1904        mpq_init(res);
1905        mpq_set_str(res,"1",10);
1906        mpq_canonicalize(res);
1907                       
1908        mpq_set_num(qkgV,kgV);
1909                       
1910//                      mpq_canonicalize(qkgV);
1911        for (int ii=0;ii<(M->colsize)-1;ii++)
1912        {
1913                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
1914                n[ii]=(int)mpz_get_d(mpq_numref(res));
1915        }
1916                        //mpz_clear(denom[this->numVars]);
1917        mpz_clear(kgV);
1918        mpq_clear(qkgV); mpq_clear(res);                       
1919                       
1920}
1921/**
1922 * For all codim-2-facets we compute the gcd of the components of the facet normal and
1923 * divide it out. Thus we get a normalized representation of each
1924 * (codim-2)-facet normal, i.e. a primitive vector
1925 */
1926void gcone::normalize()
1927{
1928        int ggT=1;
1929        facet *fAct;
1930        facet *codim2Act;
1931        fAct = this->facetPtr;
1932        codim2Act = fAct->codim2Ptr;
1933        intvec *n = new intvec(this->numVars);
1934                       
1935                        //while(codim2Act->next!=NULL)
1936        while(fAct!=NULL)
1937        {
1938                while(codim2Act!=NULL)
1939                {                               
1940                        n=codim2Act->getFacetNormal();
1941                        for(int ii=0;ii<this->numVars;ii++)
1942                        {
1943                                ggT = intgcd(ggT,(*n)[ii]);
1944                        }
1945                        for(int ii=0;ii<this->numVars;ii++)
1946                        {
1947                                (*n)[ii] = ((*n)[ii])/ggT;
1948                        }
1949                        codim2Act->setFacetNormal(n);
1950                        codim2Act = codim2Act->next;                           
1951                }
1952                fAct = fAct->next;
1953        }
1954        delete n;
1955                               
1956}
1957/** \brief Enqueue new facets into the searchlist
1958 * The searchlist (SLA for short) is implemented as a FIFO
1959 * Checks are done as follows:
1960 * 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
1961 * 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
1962 *      facet from SLA and do not add fAct.
1963 * It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
1964 * can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA.
1965 * Takes ptr to search list root
1966 * Returns a pointer to new first element of Searchlist
1967 */
1968                //void enqueueNewFacets(facet &f)
1969facet * gcone::enqueueNewFacets(facet &f)
1970{
1971        facet *slHead;
1972        slHead = &f;
1973        facet *slAct;   //called with f=SearchListRoot
1974        slAct = &f;
1975        facet *slEnd;   //Pointer to end of SLA
1976        slEnd = &f;
1977        facet *slEndStatic;     //marks the end before a new facet is added             
1978        facet *fAct;
1979        fAct = this->facetPtr;
1980        facet *codim2Act;
1981        codim2Act = this->facetPtr->codim2Ptr;
1982        facet *sl2Act;
1983        sl2Act = f.codim2Ptr;
1984        /** Pointer to a facet that will be deleted */
1985        facet *deleteMarker;
1986        deleteMarker = NULL;
1987                       
1988        /** Flag to indicate a facet that should be added to SLA*/
1989        bool doNotAdd=FALSE;
1990        /** \brief  Flag to mark a facet that might be added
1991         * The following case may occur:
1992         * A facet fAct is found to be parallel but not equal to the current facet slAct from SLA.
1993         * This does however not mean that there does not exist a facet behind the current slAct that is actually equal
1994         * to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
1995         * FALSE and the according slAct is deleted.
1996         * If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
1997         */
1998        volatile bool maybe=FALSE;
1999        /**A facet was removed, lengthOfSearchlist-- thus we must not rely on
2000         * if(notParallelCtr==lengthOfSearchList) but rather
2001         * if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
2002         */
2003        volatile bool removalOccured=FALSE;
2004        int ctr=0;      //encountered equalities in SLA
2005        int notParallelCtr=0;
2006        int lengthOfSearchList=1;
2007        while(slEnd->next!=NULL)
2008        {
2009                slEnd=slEnd->next;
2010                lengthOfSearchList++;
2011        }
2012        slEndStatic = slEnd;
2013        /*1st step: compare facetNormals*/
2014        intvec *fNormal=NULL;
2015        intvec *f2Normal=NULL;
2016        intvec *slNormal=NULL;
2017        intvec *sl2Normal=NULL;
2018                       
2019        while(fAct!=NULL)
2020        {                                               
2021                if(fAct->isFlippable==TRUE)
2022                {
2023                        maybe=FALSE;
2024                        doNotAdd=TRUE;
2025                        fNormal=fAct->getFacetNormal();
2026                        slAct = slHead;
2027                        notParallelCtr=0;
2028//                      delete deleteMarker;
2029//                      deleteMarker=NULL;
2030                        /*If slAct==NULL and fAct!=NULL
2031                        we just copy all remaining facets into SLA*/
2032                        if(slAct==NULL)
2033                        {
2034                                facet *fCopy;
2035                                fCopy = fAct;
2036                                while(fCopy!=NULL)
2037                                {
2038                                        if(slAct==NULL)
2039                                        {
2040                                                slAct = new facet();
2041                                                slHead = slAct;                                                         
2042                                        }
2043                                        else
2044                                        {
2045                                                slAct->next = new facet();
2046                                                slAct = slAct->next;
2047                                        }                                                       
2048                                        slAct->setFacetNormal(fAct->getFacetNormal());
2049                                        slAct->setUCN(fAct->getUCN());
2050                                        slAct->isFlippable=fAct->isFlippable;
2051                                        facet *f2Copy;
2052                                        f2Copy = fCopy->codim2Ptr;
2053                                        sl2Act = slAct->codim2Ptr;
2054                                        while(f2Copy!=NULL)
2055                                        {
2056                                                if(sl2Act==NULL)
2057                                                {
2058                                                        sl2Act = new facet(2);
2059                                                        slAct->codim2Ptr = sl2Act;                                     
2060                                                }
2061                                                else
2062                                                {
2063                                                        facet *marker;
2064                                                        marker = sl2Act;
2065                                                        sl2Act->next = new facet(2);
2066                                                        sl2Act = sl2Act->next;
2067                                                        sl2Act->prev = marker;
2068                                                }
2069                                                sl2Act->setFacetNormal(f2Copy->getFacetNormal());
2070                                                sl2Act->setUCN(f2Copy->getUCN());
2071                                                f2Copy = f2Copy->next;
2072                                        }
2073                                        fCopy = fCopy->next;
2074                                }
2075                                break;
2076                        }
2077                        /*End of dumping into SLA*/
2078                        while(slAct!=NULL)
2079                                        //while(slAct!=slEndStatic->next)
2080                        {
2081                                if(deleteMarker!=NULL)
2082                                {
2083//                                      delete deleteMarker;
2084//                                      deleteMarker=NULL;
2085                                }
2086                                removalOccured=FALSE;
2087                                slNormal = slAct->getFacetNormal();
2088#ifdef gfan_DEBUG
2089                                cout << "Checking facet (";
2090                                fNormal->show(1,1);
2091                                cout << ") against (";
2092                                slNormal->show(1,1);
2093                                cout << ")" << endl;
2094#endif
2095                                if(!isParallel(fNormal, slNormal))
2096                                {
2097                                        notParallelCtr++;
2098                                }
2099                                else    //fN and slN are parallel
2100                                {
2101                                        //We check whether the codim-2-facets coincide
2102                                        codim2Act = fAct->codim2Ptr;
2103                                        ctr=0;
2104                                        while(codim2Act!=NULL)
2105                                        {
2106                                                f2Normal = codim2Act->getFacetNormal();
2107                                                                //sl2Act = f.codim2Ptr;
2108                                                sl2Act = slAct->codim2Ptr;
2109                                                while(sl2Act!=NULL)
2110                                                {
2111                                                        sl2Normal = sl2Act->getFacetNormal();
2112                                                        if(areEqual(f2Normal, sl2Normal))
2113                                                                ctr++;
2114                                                        sl2Act = sl2Act->next;
2115                                                }//while(sl2Act!=NULL)
2116                                                codim2Act = codim2Act->next;
2117                                        }//while(codim2Act!=NULL)
2118                                        if(ctr==fAct->numCodim2Facets)  //facets ARE equal
2119                                        {
2120#ifdef gfan_DEBUG
2121                                                cout << "Removing facet (";
2122                                                slNormal->show(1,0);
2123                                                cout << ") from SLA:" << endl;
2124                                                fAct->fDebugPrint();
2125                                                slAct->fDebugPrint();
2126#endif
2127                                                removalOccured=TRUE;
2128                                                slAct->isFlippable=FALSE;
2129                                                doNotAdd=TRUE;
2130                                                maybe=FALSE;                                                           
2131                                                if(slAct==slHead)       //We want to delete the first element of SearchList
2132                                                {
2133                                                        deleteMarker = slAct;                           
2134                                                        slHead = slAct->next;                                           
2135                                                        if(slHead!=NULL)
2136                                                                slHead->prev = NULL;
2137//                                                      delete deleteMarker; deleteMarker=NULL;
2138                                                        //set a bool flag to mark slAct as to be deleted
2139                                                }//NOTE find a way to delete without affecting slAct = slAct->next
2140                                                else if(slAct==slEndStatic)
2141                                                {
2142                                                        deleteMarker = slAct;                                   
2143                                                        if(slEndStatic->next==NULL)
2144                                                        {                                                       
2145                                                                slEndStatic = slEndStatic->prev;
2146                                                                slEndStatic->next = NULL;
2147                                                                slEnd = slEndStatic;
2148                                                        }
2149                                                        else    //we already added a facet after slEndStatic
2150                                                        {
2151                                                                slEndStatic->prev->next = slEndStatic->next;
2152                                                                slEndStatic->next->prev = slEndStatic->prev;
2153                                                                slEndStatic = slEndStatic->prev;
2154                                                                //slEnd = slEndStatic;
2155                                                        }
2156                                                }                                                               
2157                                                else
2158                                                {
2159                                                        deleteMarker = slAct;
2160                                                        slAct->prev->next = slAct->next;
2161                                                        slAct->next->prev = slAct->prev;
2162                                                }                                                               
2163                                                //update lengthOfSearchList                                     
2164                                                lengthOfSearchList--;                                           
2165                                                break;
2166                                        }//if(ctr==fAct->numCodim2Facets)
2167                                        else    //facets are parallel but NOT equal. But this does not imply that there
2168                                                //is no other facet later in SLA that might be equal.
2169                                        {
2170                                                maybe=TRUE;//                                                                   
2171                                        }                                                       
2172                                }//if(!isParallel(fNormal, slNormal))
2173                                if( (slAct->next==NULL) && (maybe==TRUE) )
2174                                {
2175                                        doNotAdd=FALSE;
2176                                        break;
2177                                }
2178                                slAct = slAct->next;
2179                                /* NOTE The following lines must not be here but rather called inside the if-block above.
2180                                Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now,
2181                                (not nice!) but since they are in seperate branches of the if-statement there should not
2182                                be a way it gets called twice thus ommiting one facet:
2183                                slAct = slAct->next;*/                                         
2184//                                              delete deleteMarker;
2185//                                              deleteMarker=NULL;
2186                                                //if slAct was marked as to be deleted, delete it here!
2187                        }//while(slAct!=NULL)                                                                   
2188                        if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) || (doNotAdd==FALSE) )
2189                                        //if( (notParallelCtr==lengthOfSearchList ) || doNotAdd==FALSE )
2190                        {
2191#ifdef gfan_DEBUG
2192                                cout << "Adding facet (";
2193                                fNormal->show(1,0);
2194                                cout << ") to SLA " << endl;
2195#endif
2196                                                //Add fAct to SLA
2197                                facet *marker;
2198                                marker = slEnd;
2199                                facet *f2Act;
2200                                f2Act = fAct->codim2Ptr;
2201                                               
2202                                slEnd->next = new facet();
2203                                slEnd = slEnd->next;
2204                                facet *slEndCodim2Root;
2205                                facet *slEndCodim2Act;
2206                                slEndCodim2Root = slEnd->codim2Ptr;
2207                                slEndCodim2Act = slEndCodim2Root;
2208                                                               
2209                                slEnd->setUCN(this->getUCN());
2210                                slEnd->isFlippable = TRUE;
2211                                slEnd->setFacetNormal(fNormal);
2212                                slEnd->prev = marker;
2213                                //Copy codim2-facets
2214                                intvec *f2Normal;
2215                                while(f2Act!=NULL)
2216                                {
2217                                        f2Normal=f2Act->getFacetNormal();
2218                                        if(slEndCodim2Root==NULL)
2219                                        {
2220                                                slEndCodim2Root = new facet(2);
2221                                                slEnd->codim2Ptr = slEndCodim2Root;                     
2222                                                slEndCodim2Root->setFacetNormal(f2Normal);
2223                                                slEndCodim2Act = slEndCodim2Root;
2224                                        }
2225                                        else                                   
2226                                        {
2227                                                slEndCodim2Act->next = new facet(2);
2228                                                slEndCodim2Act = slEndCodim2Act->next;
2229                                                slEndCodim2Act->setFacetNormal(f2Normal);
2230                                        }
2231                                        f2Act = f2Act->next;
2232                                }
2233                                lengthOfSearchList++;                           
2234                        }//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) ||
2235                        fAct = fAct->next;
2236                }//if(fAct->isFlippable==TRUE)
2237                else
2238                {
2239                        fAct = fAct->next;
2240                }
2241        }//while(fAct!=NULL)                                           
2242        return slHead;
2243}//addC2N
2244               
2245/** \brief Compute the gcd of two ints
2246 */
2247int gcone::intgcd(int a, int b)
2248{
2249        int r, p=a, q=b;
2250        if(p < 0)
2251                p = -p;
2252        if(q < 0)
2253                q = -q;
2254        while(q != 0)
2255        {
2256                r = p % q;
2257                p = q;
2258                q = r;
2259        }
2260        return p;
2261}               
2262               
2263/** \brief Construct a dd_MatrixPtr from a cone's list of facets
2264 *
2265 */
2266dd_MatrixPtr gcone::facets2Matrix(gcone const &gc)
2267{
2268        facet *fAct;
2269        fAct = gc.facetPtr;
2270       
2271        dd_MatrixPtr M;
2272        dd_rowrange ddrows;
2273        dd_colrange ddcols;
2274        ddcols=(this->numVars)+1;
2275        ddrows=this->numFacets;
2276        dd_NumberType numb = dd_Integer;
2277        M=dd_CreateMatrix(ddrows,ddcols);                       
2278                       
2279        int jj=0;
2280       
2281        while(fAct!=NULL)
2282        {
2283                intvec *comp;
2284                comp = fAct->getFacetNormal();
2285                for(int ii=0;ii<this->numVars;ii++)
2286                {
2287                        dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
2288                }
2289                jj++;
2290                fAct=fAct->next;                               
2291        }                       
2292                       
2293        return M;
2294}
2295               
2296/** \brief Write information about a cone into a file on disk
2297 *
2298 * This methods writes the information needed for the "second" method into a file.
2299 * The file's is divided in sections containing information on
2300 * 1) the ring
2301 * 2) the cone's Groebner Basis
2302 * 3) the cone's facets
2303 * Each line contains exactly one date
2304 * Each section starts with its name in CAPITALS
2305 */
2306void gcone::writeConeToFile(gcone const &gc, bool usingIntPoints)
2307{
2308        int UCN=gc.UCN;
2309        stringstream ss;
2310        ss << UCN;
2311        string UCNstr = ss.str();               
2312                       
2313        string prefix="/tmp/cone";
2314//      string prefix="./";     //crude hack -> run tests in separate directories
2315        string suffix=".sg";
2316        string filename;
2317        filename.append(prefix);
2318        filename.append(UCNstr);
2319        filename.append(suffix);
2320       
2321       
2322//      int thisPID = getpid();
2323//      ss << thisPID;
2324//      string strPID = ss.str();
2325//      string prefix="./";
2326                                               
2327        ofstream gcOutputFile(filename.c_str());
2328        facet *fAct;
2329        fAct = gc.facetPtr;                     
2330        facet *f2Act;
2331        f2Act=fAct->codim2Ptr;
2332       
2333        char *ringString = rString(currRing);
2334                       
2335        if (!gcOutputFile)
2336        {
2337                cout << "Error opening file for writing in writeConeToFile" << endl;
2338        }
2339        else
2340        {       
2341                gcOutputFile << "UCN" << endl;
2342                gcOutputFile << this->UCN << endl;
2343                gcOutputFile << "RING" << endl; 
2344                gcOutputFile << ringString << endl;
2345                gcOutputFile << "GCBASISLENGTH" << endl;
2346                gcOutputFile << IDELEMS(gc.gcBasis) << endl;                   
2347                //Write this->gcBasis into file
2348                gcOutputFile << "GCBASIS" << endl;                             
2349                for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
2350                {                                       
2351                        gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
2352                }                               
2353                                       
2354                gcOutputFile << "FACETS" << endl;                                                               
2355               
2356                while(fAct!=NULL)
2357                {       
2358                        intvec *iv;
2359                        intvec *iv2;
2360                        iv=fAct->getFacetNormal();
2361                        f2Act=fAct->codim2Ptr;
2362                        for (int ii=0;ii<iv->length();ii++)
2363                        {
2364                                if (ii<iv->length()-1)
2365                                {
2366                                        gcOutputFile << (*iv)[ii] << ",";
2367                                }
2368                                else
2369                                {
2370                                        gcOutputFile << (*iv)[ii] << " ";
2371                                }
2372                        }
2373                        while(f2Act!=NULL)
2374                        {
2375                                iv2=f2Act->getFacetNormal();   
2376                                for(int jj=0;jj<iv2->length();jj++)
2377                                {
2378                                        if (jj<iv2->length()-1)
2379                                        {
2380                                                gcOutputFile << (*iv2)[jj] << ",";
2381                                        }
2382                                        else
2383                                        {
2384                                                gcOutputFile << (*iv2)[jj] << " ";
2385                                        }
2386                                }
2387                                f2Act = f2Act->next;
2388                        }
2389                        gcOutputFile << endl;
2390                        fAct=fAct->next;
2391                }                       
2392        gcOutputFile.close();
2393        }
2394                       
2395}//writeConeToFile(gcone const &gc)
2396               
2397/** \brief Reads a cone from a file identified by its number
2398 */
2399void gcone::readConeFromFile(int UCN, gcone *gc)
2400{
2401        //int UCN=gc.UCN;
2402        stringstream ss;
2403        ss << UCN;
2404        string UCNstr = ss.str();
2405        string line;
2406        string strGcBasisLength;
2407        string strMonom, strCoeff, strCoeffNom, strCoeffDenom;         
2408        int gcBasisLength=0;
2409        int intCoeff=1;
2410        int intCoeffNom=1;              //better (number) but that's not compatible with stringstream;
2411        int intCoeffDenom=1;
2412        number nCoeff=nInit(1);
2413        number nCoeffNom=nInit(1);
2414        number nCoeffDenom=nInit(1);
2415        bool hasCoeffInQ = FALSE;       //does the polynomial have rational coeff?
2416        bool hasNegCoeff = FALSE;       //or a negative one?
2417        size_t found;                   //used for find_first_(not)_of
2418
2419        string prefix="/tmp/cone";
2420        string suffix=".sg";
2421        string filename;
2422        filename.append(prefix);
2423        filename.append(UCNstr);
2424        filename.append(suffix);
2425                                       
2426        ifstream gcInputFile(filename.c_str(), ifstream::in);
2427       
2428        ring saveRing=currRing; 
2429        rChangeCurrRing(gc->baseRing);
2430        poly strPoly=pInit();
2431        poly resPoly=pInit();   //The poly to be read in
2432       
2433        string::iterator EOL;
2434        int terms=1;    //#Terms in the poly
2435       
2436        while( !gcInputFile.eof() )
2437        {
2438                getline(gcInputFile,line);
2439                hasCoeffInQ = FALSE;
2440                hasNegCoeff = FALSE;
2441               
2442                if(line=="GCBASISLENGTH")
2443                {
2444                        getline(gcInputFile, line);
2445                        strGcBasisLength = line;
2446                        gcBasisLength=atoi(strGcBasisLength.c_str());
2447                }
2448                if(line=="GCBASIS")
2449                {
2450                        for(int jj=0;jj<gcBasisLength;jj++)
2451                        {
2452                                getline(gcInputFile,line);
2453                                //find out how many terms the poly consists of
2454                                for(EOL=line.begin(); EOL!=line.end();++EOL)
2455                                {
2456                                        string s;
2457                                        s=*EOL;
2458                                        if(s=="+" || s=="-")
2459                                                terms++;
2460                                }
2461                                //magically convert strings into polynomials
2462                                //polys.cc:p_Read
2463                                //check until first occurance of + or -
2464                                //data or c_str
2465                                while(!line.empty())
2466                                {
2467//                                      resPoly=pInit();
2468//                                      strPoly=pInit();
2469                                        hasNegCoeff = FALSE;
2470                                        hasCoeffInQ = FALSE;
2471                                        found = line.find_first_of("+-");       //get the first monomial
2472                                        string tmp;
2473                                        tmp=line[found];
2474//                                      if(found!=0 && (tmp.compare("-")==0) )
2475//                                              hasNegCoeff = TRUE;     //We have a coeff < 0
2476                                        if(found==0)
2477                                        {
2478                                                if(tmp.compare("-")==0)
2479                                                        hasNegCoeff = TRUE;
2480                                                line.erase(0,1);        //remove leading + or -
2481                                                found = line.find_first_of("+-");       //adjust found
2482                                        }
2483                                        strMonom = line.substr(0,found);
2484                                        line.erase(0,found);
2485                                        found = strMonom.find_first_of("/");
2486                                        if(found!=string::npos) //i.e. "/" exists in strMonom
2487                                        {
2488                                                hasCoeffInQ = TRUE;
2489                                                strCoeffNom=strMonom.substr(0,found);                                           
2490                                                strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890",found+1));
2491                                                strMonom.erase(0,found);
2492                                                strMonom.erase(0,strMonom.find_first_not_of("1234567890/"));                   
2493//                                              ss << strCoeffNom;
2494//                                              ss >> intCoeffNom;
2495//                                              nCoeffNom=(snumber*)strCoeffNom.c_str();
2496                                                nRead(strCoeffNom.c_str(), &nCoeffNom);
2497                                                nRead(strCoeffDenom.c_str(), &nCoeffDenom);
2498//                                              nCoeffNom=nInit(intCoeffNom);
2499//                                              ss << strCoeffDenom;
2500//                                              ss >> intCoeffDenom;
2501//                                              nCoeffDenom=nInit(intCoeffDenom);
2502//                                              nCoeffDenom=(snumber*)strCoeffNom.c_str();
2503                                                //NOTE NOT SURE WHETHER THIS WILL WORK!
2504                                                //nCoeffNom=nInit(intCoeffNom);
2505//                                              nCoeffDenom=(number)strCoeffDenom.c_str();
2506//                                              nCoeffDenom=(number)strCoeffDenom.c_str();
2507                                        }
2508                                        else
2509                                        {
2510                                                found = strMonom.find_first_not_of("1234567890");
2511                                                strCoeff = strMonom.substr(0,found);
2512                                                if(!strCoeff.empty())
2513                                                {
2514                                                        nRead(strCoeff.c_str(),&nCoeff);
2515                                                }
2516                                                else
2517                                                {
2518//                                                      intCoeff = 1;
2519                                                        nCoeff = nInit(1);
2520                                                }
2521                                                                                               
2522                                        }
2523                                        const char* monom = strMonom.c_str();
2524                                               
2525                                        p_Read(monom,strPoly,currRing);
2526                                        switch (hasCoeffInQ)
2527                                        {
2528                                                case TRUE:
2529                                                        if(hasNegCoeff)
2530                                                                nCoeffNom=nNeg(nCoeffNom);
2531//                                                              intCoeffNom *= -1;
2532//                                                      pSetCoeff(strPoly, nDiv((number)intCoeffNom, (number)intCoeffDenom));
2533                                                        pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom));
2534                                                        break;
2535                                                case FALSE:
2536                                                        if(hasNegCoeff)
2537                                                                nCoeff=nNeg(nCoeff);
2538//                                                              intCoeff *= -1;
2539                                                        if(!nIsOne(nCoeff))
2540                                                        {
2541//                                                              if(hasNegCoeff)
2542//                                                                      intCoeff *= -1;
2543//                                                              pSetCoeff(strPoly,(number) intCoeff);
2544                                                                pSetCoeff(strPoly, nCoeff );
2545                                                        }
2546                                                        break;
2547                                                                                                       
2548                                        }
2549                                                //pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1???
2550                                        if(pIsConstantComp(resPoly))
2551                                                resPoly=pCopy(strPoly);                                                 
2552                                        else
2553                                                resPoly=pAdd(resPoly,strPoly);
2554                                       
2555                                       
2556                                }//while(!line.empty())         
2557                       
2558                                gc->gcBasis->m[jj]=pCopy(resPoly);
2559                                resPoly=NULL;   //reset
2560                        }
2561                        break;
2562                }               
2563        }//while(!gcInputFile.eof())
2564       
2565        gcInputFile.close();
2566        rChangeCurrRing(saveRing);
2567}
2568       
2569// static void gcone::idPrint(ideal &I)
2570// {
2571//      for(int ii=0;ii<IDELEMS(I);ii++)
2572//      {
2573//              cout << "_[" << ii << "]=" << I->m[ii] << endl;
2574//      }
2575// }
2576ring gcone::getBaseRing()
2577{
2578        return rCopy(this->baseRing);
2579}
2580/** \brief Gather the output
2581* List of lists
2582*\param *gc Pointer to gcone, preferably gcRoot ;-)
2583*\param n the number of cones
2584*
2585*/
2586lists lprepareResult(gcone *gc, int n)
2587{
2588        gcone *gcAct;
2589        gcAct = gc;     
2590        facet *fAct;
2591        fAct = gc->facetPtr;
2592       
2593        lists res=(lists)omAllocBin(slists_bin);
2594        res->Init(n);   //initialize to store n cones
2595        for(int ii=0;ii<n;ii++)
2596        {
2597                res->m[ii].rtyp=LIST_CMD;
2598                lists l=(lists)omAllocBin(slists_bin);
2599                l->Init(6);
2600                l->m[0].rtyp=INT_CMD;
2601                l->m[0].data=(void*)gcAct->getUCN();
2602                l->m[1].rtyp=IDEAL_CMD;
2603                /*The following is necessary for leaves in the tree of cones
2604                * Since we don't use them in the computation and gcBasis is
2605                * set to (poly)NULL in noRevS we need to get this back here.
2606                */
2607                if(gcAct->gcBasis->m[0]==(poly)NULL)
2608                        gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
2609                l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
2610
2611                l->m[2].rtyp=INTVEC_CMD;
2612                intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));
2613                l->m[2].data=(void*)ivCopy(&iv);
2614               
2615                l->m[3].rtyp=LIST_CMD;
2616                        lists lCodim2List = (lists)omAllocBin(slists_bin);
2617                        lCodim2List->Init(gcAct->numFacets);           
2618 
2619                        fAct = gcAct->facetPtr;//fAct->codim2Ptr;
2620                        int jj=0;
2621                        while(fAct!=NULL && jj<gcAct->numFacets)
2622                        {
2623                                lCodim2List->m[jj].rtyp=INTVEC_CMD;
2624                                intvec ivC2=(gcAct->f2M(gcAct,fAct,2));
2625                                lCodim2List->m[jj].data=(void*)ivCopy(&ivC2);
2626                                jj++;
2627                                fAct = fAct->next;
2628                        }               
2629                l->m[3].data=(void*)lCodim2List;
2630               
2631               
2632                l->m[4].rtyp=RING_CMD;
2633                l->m[4].data=(void*)(gcAct->getBaseRing());             
2634                l->m[5].rtyp=INT_CMD;
2635                l->m[5].data=(void*)gcAct->getPredUCN();
2636                res->m[ii].data=(void*)l;
2637                gcAct = gcAct->next;
2638        }               
2639       
2640        return res;
2641}
2642/** \brief Write facets of a cone into a matrix
2643* Takes a pointer to a facet as 2nd arg
2644* f should always point to gc->facetPtr
2645* param n is used to determine whether it operates in codim 1 or 2
2646*
2647*/
2648intvec gcone::f2M(gcone *gc, facet *f, int n)
2649{
2650        facet *fAct;
2651        intvec *res;
2652        intvec *fNormal = new intvec(gc->numVars);
2653//      int codim=n;
2654//      int bound;
2655//      if(f==gc->facetPtr)
2656        if(n==1)
2657        {
2658                intvec *m1Res=new intvec(gc->numFacets,gc->numVars,0);
2659                res = m1Res;
2660                fAct = gc->facetPtr;
2661//              bound = gc->numFacets*(this->numVars);         
2662        }
2663        else
2664        {
2665                fAct = f->codim2Ptr;
2666                intvec *m2Res = new intvec(f->numCodim2Facets,gc->numVars,0);
2667                res = m2Res;           
2668//              bound = fAct->numCodim2Facets*(this->numVars);
2669
2670        }
2671        int ii=0;
2672        while(fAct!=NULL )//&& ii < bound )
2673        {
2674                fNormal = fAct->getFacetNormal();
2675                for(int jj=0;jj<this->numVars;jj++)
2676                {
2677                        (*res)[ii]=(*fNormal)[jj];
2678                        ii++;
2679                }
2680                fAct = fAct->next;
2681        }
2682
2683        delete fNormal;
2684        return *res;
2685}
2686
2687int gcone::counter=0;
2688int gfanHeuristic;
2689// ideal gfan(ideal inputIdeal, int h)
2690lists gfan(ideal inputIdeal, int h)
2691{
2692        lists lResList; //this is the object we return
2693       
2694        if(rHasGlobalOrdering(currRing))
2695        {       
2696        int numvar = pVariables; 
2697        gfanHeuristic = h;
2698       
2699        enum searchMethod {
2700                reverseSearch,
2701                noRevS
2702        };
2703       
2704        searchMethod method;
2705        method = noRevS;
2706
2707       
2708#ifdef gfan_DEBUG
2709        cout << "Now in subroutine gfan" << endl;
2710#endif
2711        ring inputRing=currRing;        // The ring the user entered
2712        ring rootRing;                  // The ring associated to the target ordering
2713        ideal res;     
2714        facet *fRoot;   
2715       
2716        if(method==noRevS)
2717        {
2718                gcone *gcRoot = new gcone(currRing,inputIdeal);
2719                gcone *gcAct;
2720                gcAct = gcRoot;
2721                gcAct->numVars=pVariables;
2722                gcAct->getGB(inputIdeal);               
2723#ifndef NDEBUG
2724                cout << "GB of input ideal is:" << endl;
2725                idShow(gcAct->gcBasis);
2726#endif
2727                if(gcAct->isMonomial(gcAct->gcBasis))
2728                {
2729                        WerrorS("Monomial input - terminating");
2730                        res=gcAct->gcBasis;
2731                        //break;
2732                }
2733                cout << endl;
2734                gcAct->getConeNormals(gcAct->gcBasis);         
2735                gcAct->noRevS(*gcAct);         
2736                //res=gcAct->gcBasis;
2737                //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
2738                res = inputIdeal; 
2739                lResList=lprepareResult(gcRoot,gcRoot->getCounter());
2740                /*Cleanup*/
2741                gcone *gcDel;   
2742                gcDel = gcRoot;
2743                gcAct = gcRoot;
2744                while(gcAct!=NULL)
2745                {
2746                        gcDel = gcAct;
2747                        gcAct = gcAct->next;
2748                        delete gcDel;
2749                }
2750        }
2751        dd_free_global_constants();
2752        }//rHasGlobalOrdering
2753        else
2754        {
2755                WerrorS("Ring has non-global ordering - terminating");
2756                lResList->Init(1);
2757                lResList->m[0].rtyp=INT_CMD;
2758                int ires=0;
2759                lResList->m[0].data=(void*)&ires;
2760        }
2761       
2762        /*Return result*/
2763        return lResList;
2764}
2765
2766#endif
Note: See TracBrowser for help on using the repository browser.