source: git/kernel/gfan.cc @ 2e658f

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