source: git/kernel/gfan.cc @ 5684df

spielwiese
Last change on this file since 5684df was 5684df, checked in by Martin Monerjan, 15 years ago
found bug in flippability test git-svn-id: file:///usr/local/Singular/svn/trunk@11767 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 38.3 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009-05-04 15:14:38 $
5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.44 2009-05-04 15:14:38 monerjan Exp $
6$Id: gfan.cc,v 1.44 2009-05-04 15:14:38 monerjan Exp $
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 "prCopy.h"
23#include <iostream.h>   //deprecated
24#include <bitset>
25
26/*remove the following at your own risk*/
27#ifndef GMPRATIONAL
28#define GMPRATIONAL
29#endif
30
31//Hacks for different working places
32#define ITWM
33
34#ifdef UNI
35#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
36#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
37#endif
38
39#ifdef HOME
40#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
41#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
42#endif
43
44#ifdef ITWM
45#include "/u/slg/monerjan/cddlib/include/setoper.h"
46#include "/u/slg/monerjan/cddlib/include/cdd.h"
47#include "/u/slg/monerjan/cddlib/include/cddmp.h"
48#endif
49
50#ifndef gfan_DEBUG
51#define gfan_DEBUG
52#endif
53
54//#include gcone.h
55
56/**
57*\brief Class facet
58*       Implements the facet structure as a linked list
59*
60*/
61class facet
62{
63        private:
64                /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
65                intvec *fNormal;
66               
67                /** \brief The Groebner basis on the other side of a shared facet
68                 *
69                 * In order not to have to compute the flipped GB twice we store the basis we already get
70                 * when identifying search facets. Thus in the next step of the reverse search we can
71                 * just copy the old cone and update the facet and the gcBasis.
72                 * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
73                 */
74                ideal flipGB;           //The Groebner Basis on the other side, computed via gcone::flip
75
76                       
77        public:
78                /** The default constructor. Do I need a constructor of type facet(intvec)? */
79                facet()                 
80                {
81                        // Pointer to next facet.  */
82                        /* Defaults to NULL. This way there is no need to check explicitly */
83                        this->next=NULL; 
84                }
85               
86                /** The default destructor */
87                ~facet(){;}
88               
89                /** Stores the facet normal \param intvec*/
90                void setFacetNormal(intvec *iv)
91                {
92                        this->fNormal = ivCopy(iv);
93                        //return;
94                }
95               
96                /** Hopefully returns the facet normal */
97                intvec *getFacetNormal()
98                {       
99                        //this->fNormal->show();        cout << endl;   
100                        return this->fNormal;
101                }
102               
103                /** Method to print the facet normal*/
104                void printNormal()
105                {
106                        fNormal->show();
107                }
108               
109                /** Store the flipped GB*/
110                void setFlipGB(ideal I)
111                {
112                        this->flipGB=I;
113                }
114               
115                /** Return the flipped GB*/
116                ideal getFlipGB()
117                {
118                        return this->flipGB;
119                }
120               
121                /** Print the flipped GB*/
122                void printFlipGB()
123                {
124                        idShow(this->flipGB);
125                }
126               
127                bool isFlippable;       //flippable facet? Want to have cone->isflippable.facet[i]
128                bool isIncoming;        //Is the facet incoming or outgoing?
129                facet *next;            //Pointer to next facet
130               
131                friend class gcone;     //Bad style
132};
133
134/**
135*\brief Implements the cone structure
136*
137* A cone is represented by a linked list of facet normals
138* @see facet
139*/
140/*class gcone
141finally this should become s.th. like gconelib.{h,cc} to provide an API
142*/
143class gcone
144{
145        private:
146                int numFacets;          //#of facets of the cone
147                ring rootRing;          //good to know this -> generic walk
148                ideal inputIdeal;       //the original
149                ring baseRing;          //the basering of the cone             
150                /* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
151                intvec *ivIntPt;        //an interior point of the cone
152               
153        public:
154                /** \brief Default constructor.
155                *
156                * Initialises this->next=NULL and this->facetPtr=NULL
157                */
158                gcone()
159                {
160                        this->next=NULL;
161                        this->facetPtr=NULL;
162                        this->baseRing=currRing;
163                }
164               
165                /** \brief Constructor with ring and ideal
166                *
167                * This constructor takes the root ring and the root ideal as parameters and stores
168                * them in the private members gcone::rootRing and gcone::inputIdeal
169                */
170                gcone(ring r, ideal I)
171                {
172                        this->next=NULL;
173                        this->facetPtr=NULL;
174                        this->rootRing=r;
175                        this->inputIdeal=I;
176                        this->baseRing=currRing;
177                }
178               
179                /** \brief Copy constructor
180                *
181                * Copies one cone, sets this->gcBasis to the flipped GB and reverses the
182                * direction of the according facet normal
183                */                     
184                gcone(const gcone& gc)         
185                {
186                        this->next=NULL;
187                        this->numVars=gc.numVars;                       
188                        facet *fAct= new facet();                       
189                        this->facetPtr=fAct;
190                       
191                        intvec *ivtmp = new intvec(this->numVars);                                             
192                        ivtmp = gc.facetPtr->getFacetNormal();
193                        ivtmp->show();                 
194                       
195                        ideal gb;
196                        gb=gc.facetPtr->getFlipGB();                   
197                        this->gcBasis=gb;//gc.facetPtr->getFlipGB();    //this cone's GB is the flipped GB
198                        idShow(this->gcBasis);
199                       
200                        /*Reverse direction of the facet normal to make it an inner normal*/                   
201                        for (int ii=0; ii<this->numVars;ii++)
202                        {
203                                (*ivtmp)[ii]=-(*ivtmp)[ii];                             
204                        }
205                        //ivtmp->show(); cout << endl;
206                        fAct->setFacetNormal(ivtmp);
207                        //fAct->printNormal();cout << endl;
208                        //ivtmp->show();cout << endl;
209                }
210               
211                /** \brief Default destructor */
212                ~gcone(){;}             //destructor
213               
214                /** Pointer to the first facet */
215                facet *facetPtr;        //Will hold the adress of the first facet; set by gcone::getConeNormals
216               
217                /** # of variables in the ring */
218                int numVars;            //#of variables in the ring
219               
220                /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
221                ideal gcBasis;          //GB of the cone, set by gcone::getGB();
222                gcone *next;            //Pointer to *previous* cone in search tree
223               
224                /** \brief Set the interior point of a cone */
225                void setIntPoint(intvec *iv)
226                {
227                        this->ivIntPt=ivCopy(iv);
228                }
229               
230                /** \brief Return the interior point */
231                intvec *getIntPoint()
232                {
233                        return this->ivIntPt;
234                }
235               
236                void showIntPoint()
237                {
238                        ivIntPt->show();
239                }
240               
241                /** \brief Compute the normals of the cone
242                *
243                * This method computes a representation of the cone in terms of facet normals. It takes an ideal
244                * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
245                * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
246                * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
247                * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
248                * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
249                * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
250                *
251                * Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
252                * an interior point of the cone.
253                */
254                void getConeNormals(ideal const &I, bool compIntPoint=FALSE)
255                {
256#ifdef gfan_DEBUG
257                        std::cout << "*** Computing Inequalities... ***" << std::endl;
258#endif         
259                        //All variables go here - except ineq matrix and *v, see below
260                        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
261                        int pCompCount;                 // # of terms in a poly
262                        poly aktpoly;
263                        int numvar = pVariables;        // # of variables in a polynomial (or ring?)
264                        int leadexp[numvar];            // dirty hack of exp.vects
265                        int aktexp[numvar];
266                        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
267                        dd_rowrange ddrows;
268                        dd_colrange ddcols;
269                        dd_rowset ddredrows;            // # of redundant rows in ddineq
270                        dd_rowset ddlinset;             // the opposite
271                        dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
272                        dd_NumberType ddnumb=dd_Integer; //Number type
273                        dd_ErrorType dderr=dd_NoError;  //
274                        // End of var declaration
275#ifdef gfan_DEBUG
276                        cout << "The Groebner basis has " << lengthGB << " elements" << endl;
277                        cout << "The current ring has " << numvar << " variables" << endl;
278#endif
279                        cols = numvar;
280               
281                        //Compute the # inequalities i.e. rows of the matrix
282                        rows=0; //Initialization
283                        for (int ii=0;ii<IDELEMS(I);ii++)
284                        {
285                                aktpoly=(poly)I->m[ii];
286                                rows=rows+pLength(aktpoly)-1;
287                        }
288#ifdef gfan_DEBUG
289                        cout << "rows=" << rows << endl;
290                        cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
291#endif
292                        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
293                        dd_set_global_constants();
294                        ddrows=rows;
295                        ddcols=cols;
296                        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
297                        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
298               
299                        // We loop through each g\in GB and compute the resulting inequalities
300                        for (int i=0; i<IDELEMS(I); i++)
301                        {
302                                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
303                                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
304                                cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
305               
306                                int *v=(int *)omAlloc((numvar+1)*sizeof(int));
307                                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
308                               
309                                //Store leadexp for aktpoly
310                                for (int kk=0;kk<numvar;kk++)
311                                {
312                                        leadexp[kk]=v[kk+1];
313                                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
314                                        //but compute the diff below
315                                }
316               
317                               
318                                while (pNext(aktpoly)!=NULL) //move to next term until NULL
319                                {
320                                        aktpoly=pNext(aktpoly);
321                                        pSetm(aktpoly);         //doesn't seem to help anything
322                                        pGetExpV(aktpoly,v);
323                                        for (int kk=0;kk<numvar;kk++)
324                                        {
325                                                aktexp[kk]=v[kk+1];
326                                                //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
327                                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
328                                        }
329                                        aktmatrixrow=aktmatrixrow+1;
330                                }//while
331               
332                        } //for
333               
334                        //Maybe add another row to contain the constraints of the standard simplex?
335
336#ifdef gfan_DEBUG
337                        cout << "The inequality matrix is" << endl;
338                        dd_WriteMatrix(stdout, ddineq);
339#endif
340
341                        // The inequalities are now stored in ddineq
342                        // Next we check for superflous rows
343                        ddredrows = dd_RedundantRows(ddineq, &dderr);
344                        if (dderr!=dd_NoError)                  // did an error occur?
345                        {
346                                dd_WriteErrorMessages(stderr,dderr);    //if so tell us
347                        } else
348                        {
349                                cout << "Redundant rows: ";
350                                set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
351                        }//if dd_Error
352               
353                        //Remove reduntant rows here!
354                        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
355                        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
356                        ddcols = ddineq->colsize;
357#ifdef gfan_DEBUG
358                        cout << "Having removed redundancies, the normals now read:" << endl;
359                        dd_WriteMatrix(stdout,ddineq);
360                        cout << "Rows = " << ddrows << endl;
361                        cout << "Cols = " << ddcols << endl;
362#endif
363                       
364                        /*Write the normals into class facet*/
365#ifdef gfan_DEBUG
366                        cout << "Creating list of normals" << endl;
367#endif
368                        /*The pointer *fRoot should be the return value of this function*/
369                        facet *fRoot = new facet();     //instantiate new facet
370                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
371                        facet *fAct;                    //instantiate pointer to active facet
372                        fAct = fRoot;                   //Seems to do the trick. fRoot and fAct have to point to the same adress!
373                        std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
374                        for (int kk = 0; kk<ddrows; kk++)
375                        {
376                                intvec *load = new intvec(numvar);      //intvec to store a single facet normal that will then be stored via setFacetNormal
377                                for (int jj = 1; jj <ddcols; jj++)
378                                {
379#ifdef GMPRATIONAL
380                                        double foo;
381                                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
382/*#ifdef gfan_DEBUG
383                                        std::cout << "fAct is " << foo << " at " << fAct << std::endl;
384#endif*/
385                                        (*load)[jj-1] = (int)foo;               //store typecasted entry at pos jj-1 of load
386#endif
387#ifndef GMPRATIONAL
388                                        double *foo;
389                                        foo = (double*)ddineq->matrix[kk][jj];  //get entry from actual position#endif
390/*#ifdef gfan_DEBUG
391                                        std::cout << "fAct is " << *foo << " at " << fAct << std::endl;
392#endif*/
393                                        (*load)[jj-1] = (int)*foo;              //store typecasted entry at pos jj-1 of load
394#endif //GMPRATIONAL                                   
395                       
396
397                                        //(*load)[jj-1] = (int)foo;             //store typecasted entry at pos jj-1 of load
398                                        //check for flipability here
399                                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
400                                        {
401                                                //fAct->next = new facet();     //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
402                                        }
403                                }//for (int jj = 1; jj <ddcols; jj++)
404                                /*Quick'n'dirty hack for flippability*/ 
405                                bool isFlippable;
406                                //NOTE BUG HERE
407                                /* The flippability check is wrong!
408                                (1,-4) will pass, but (-1,7) will not.
409                                REWRITE THIS
410                                */                     
411                                for (int jj = 0; jj<this->numVars; jj++)
412                                {                                       
413                                        intvec *ivCanonical = new intvec(this->numVars);
414                                        (*ivCanonical)[jj]=1;                                   
415                                        if (dotProduct(load,ivCanonical)>=0)
416                                        {
417                                                isFlippable=FALSE;                                             
418                                        }
419                                        else
420                                        {
421                                                isFlippable=TRUE;
422                                        }                                       
423                                        delete ivCanonical;
424                                }//for (int jj = 0; jj<this->numVars; jj++)
425                                if (isFlippable==FALSE)
426                                {
427                                        cout << "Ignoring facet";
428                                        load->show();
429                                        //fAct->next=NULL;
430                                }
431                                else
432                                {       /*Now load should be full and we can call setFacetNormal*/
433                                        fAct->setFacetNormal(load);
434                                        fAct->next = new facet();
435                                        //fAct->printNormal();
436                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
437                                }//if (isFlippable==FALSE)
438                                delete load;
439                        }//for (int kk = 0; kk<ddrows; kk++)
440                        /*
441                        Now we should have a linked list containing the facet normals of those facets that are
442                        -irredundant
443                        -flipable
444                        Adressing is done via *facetPtr
445                        */
446                       
447                        if (compIntPoint==TRUE)
448                        {
449                                intvec *iv = new intvec(this->numVars);
450                                interiorPoint(ddineq, *iv);
451                                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
452                                delete iv;
453                        }
454                       
455                        //Clean up but don't delete the return value! (Whatever it will turn out to be)
456                       
457                        dd_FreeMatrix(ddineq);
458                        set_free(ddredrows);
459                        free(ddnewpos);
460                        set_free(ddlinset);
461                        dd_free_global_constants();
462
463                }//method getConeNormals(ideal I)       
464               
465               
466                /** \brief Compute the Groebner Basis on the other side of a shared facet
467                *
468                * Implements algorithm 4.3.2 from Anders' thesis.
469                * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
470                * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
471                * is parallel to \f$ leadexp(g) \f$
472                * Parallelity is checked using basic linear algebra. See gcone::isParallel.
473                * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
474                * computing an interior point of the facet and taking all terms having the same weight with respect
475                * to this interior point.
476                *\param ideal, facet
477                * Input: a marked,reduced Groebner basis and a facet
478                */
479                void flip(ideal gb, facet *f)           //Compute "the other side"
480                {                       
481                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
482                        fNormal = f->getFacetNormal();  //read this->fNormal;
483#ifdef gfan_DEBUG
484                        std::cout << "===" << std::endl;
485                        std::cout << "running gcone::flip" << std::endl;
486                        std::cout << "fNormal=";
487                        fNormal->show();
488                        std::cout << std::endl;
489#endif                         
490                        /*1st step: Compute the initial ideal*/
491                        poly initialFormElement[IDELEMS(gb)];   //array of #polys in GB to store initial form
492                        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
493                        poly aktpoly;
494                        intvec *check = new intvec(this->numVars);      //array to store the difference of LE and v
495                       
496                        for (int ii=0;ii<IDELEMS(gb);ii++)
497                        {
498                                aktpoly = (poly)gb->m[ii];                                                             
499                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
500                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
501                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
502                                initialFormElement[ii]=pHead(aktpoly);
503                               
504                                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
505                                {
506                                        aktpoly=pNext(aktpoly); //next term
507                                        pSetm(aktpoly);
508                                        pGetExpV(aktpoly,v);           
509                                        /* Convert (int)v into (intvec)check */                 
510                                        for (int jj=0;jj<this->numVars;jj++)
511                                        {
512                                                //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
513                                                //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
514                                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
515                                        }
516#ifdef gfan_DEBUG
517                                        cout << "check=";                       
518                                        check->show();
519                                        cout << endl;
520#endif
521                                        //TODO why not *check, *fNormal????
522                                        if (isParallel(*check,*fNormal)) //pass *check when
523                                        {
524                                                cout << "Parallel vector found, adding to initialFormElement" << endl;                 
525                                                initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
526                                        }                                               
527                                }//while
528#ifdef gfan_DEBUG
529                                cout << "Initial Form=";                               
530                                pWrite(initialFormElement[ii]);
531                                cout << "---" << endl;
532#endif
533                                /*Now initialFormElement must be added to (ideal)initialForm */
534                                initialForm->m[ii]=initialFormElement[ii];
535                        }//for
536                        //f->setFlipGB(initialForm);    //FIXME PROBABLY WRONG TO STORE HERE SINCE INA!=flibGB                 
537#ifdef gfan_DEBUG
538                        cout << "Initial ideal is: " << endl;
539                        idShow(initialForm);
540                        //f->printFlipGB();
541                        cout << "===" << endl;
542#endif
543                        delete check;
544                       
545                        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
546                        /*Substep 2.1
547                        compute $G_{-\alpha}(in_v(I))
548                        see journal p. 66
549                        */
550                        ring srcRing=currRing;
551                       
552                        /* copied and modified from ring.cc::rAssureSyzComp */
553                        //ring tmpRing=rCopyAndChangeWeight(srcRing,fNormal);
554                        ring tmpRing=rCopy0(srcRing);
555                        int i=rBlocks(srcRing);
556                        int j;
557                        tmpRing->order=(int *)omAlloc((i+1)*sizeof(int));
558                        /*NOTE This should probably be set, but as of now crashes Singular*/
559                        /*tmpRing->block0=(int *)omAlloc0((i+1)*sizeof(int));
560                        tmpRing->block1=(int *)omAlloc0((i+1)*sizeof(int));*/
561                        tmpRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
562                        for(j=i;j>0;j--)
563                        {
564                                tmpRing->order[j]=srcRing->order[j-1];
565                                tmpRing->block0[j]=srcRing->block0[j-1];
566                                tmpRing->block1[j]=srcRing->block1[j-1];
567                                if (srcRing->wvhdl[j-1] != NULL)
568                                {
569                                        tmpRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
570                                }
571                        }
572                        tmpRing->order[0]=ringorder_a;
573                        tmpRing->order[1]=ringorder_dp;
574                        tmpRing->order[2]=ringorder_C;
575                        //tmpRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
576                       
577                        for (int ii=0;ii<this->numVars;ii++)
578                        {                               
579                                tmpRing->wvhdl[0][ii]=-(*fNormal)[ii];  //What exactly am I doing here?
580                                //cout << tmpring->wvhdl[0][ii] << endl;
581                        }
582                        rComplete(tmpRing);
583                        rChangeCurrRing(tmpRing);
584                       
585                        rWrite(currRing); cout << endl;
586                        ideal ina;                     
587                        ina=idrCopyR(initialForm,srcRing);                     
588#ifdef gfan_DEBUG
589                        cout << "ina=";
590                        idShow(ina); cout << endl;
591#endif
592                        ideal H;
593                        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
594                        idSkipZeroes(H);
595#ifdef gfan_DEBUG
596                        cout << "H="; idShow(H); cout << endl;
597#endif
598                        /*Substep 2.2
599                        do the lifting and mark according to H
600                        */
601                        rChangeCurrRing(srcRing);
602                        ideal srcRing_H;
603                        ideal srcRing_HH;                       
604                        srcRing_H=idrCopyR(H,tmpRing);
605#ifdef gfan_DEBUG
606                        cout << "srcRing_H = ";
607                        idShow(srcRing_H); cout << endl;
608#endif
609                        srcRing_HH=ffG(srcRing_H,this->gcBasis);               
610#ifdef gfan_DEBUG
611                        cout << "srcRing_HH = ";
612                        idShow(srcRing_HH); cout << endl;
613#endif
614                        /*Substep 2.2.1
615                        Mark according to G_-\alpha
616                        Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
617                        we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
618                        represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
619                        Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
620                        compute the difference accordingly
621                        */
622                        dd_set_global_constants();
623                        bool markingsAreCorrect=FALSE;
624                        dd_MatrixPtr intPointMatrix;
625                        int iPMatrixRows=0;
626                        dd_rowrange aktrow=0;                   
627                        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
628                        {
629                                poly aktpoly=(poly)srcRing_HH->m[ii];
630                                iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
631                        }
632                        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
633                        construction of the differences
634                        */
635                        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
636                        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
637                       
638                        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
639                        {
640                                markingsAreCorrect=FALSE;       //crucial to initialise here
641                                poly aktpoly=srcRing_HH->m[ii];
642                                /*Comparison of leading monomials is done via exponent vectors*/
643                                for (int jj=0;jj<IDELEMS(H);jj++)
644                                {
645                                        int *src_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
646                                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
647                                        pGetExpV(aktpoly,src_ExpV);
648                                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
649                                        pGetExpV(pCopy(H->m[ii]),dst_ExpV);
650                                        rChangeCurrRing(srcRing);
651                                        bool expVAreEqual=TRUE;
652                                        for (int kk=1;kk<=this->numVars;kk++)
653                                        {
654                                                cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
655                                                if (src_ExpV[kk]!=dst_ExpV[kk])
656                                                {
657                                                        expVAreEqual=FALSE;
658                                                }
659                                        }                                       
660                                        //if (*src_ExpV == *dst_ExpV)
661                                        if (expVAreEqual==TRUE)
662                                        {
663                                                markingsAreCorrect=TRUE; //everything is fine
664                                                cout << "correct markings" << endl;
665                                        }//if (pHead(aktpoly)==pHead(H->m[jj])
666                                        delete src_ExpV;
667                                        delete dst_ExpV;
668                                }//for (int jj=0;jj<IDELEMS(H);jj++)
669                               
670                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
671                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
672                                if (markingsAreCorrect==TRUE)
673                                {
674                                        pGetExpV(aktpoly,leadExpV);
675                                }
676                                else
677                                {
678                                        rChangeCurrRing(tmpRing);
679                                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
680                                        rChangeCurrRing(srcRing);
681                                }
682                                /*compute differences of the expvects*/                         
683                                while (pNext(aktpoly)!=NULL)
684                                {
685                                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
686                                        is not omitted when computing the differences*/
687                                        if(markingsAreCorrect==TRUE)
688                                        {
689                                                aktpoly=pNext(aktpoly);
690                                                pGetExpV(aktpoly,v);
691                                        }
692                                        else
693                                        {
694                                                pGetExpV(pHead(aktpoly),v);
695                                                markingsAreCorrect=TRUE;
696                                        }
697
698                                        for (int jj=0;jj<this->numVars;jj++)
699                                        {                               
700                                                /*Store into ddMatrix*/                                                         
701                                                dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
702                                        }
703                                        aktrow +=1;
704                                }
705                                delete v;
706                                delete leadExpV;
707                        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
708                        /*Now we add the constraint for the standard simplex*/
709                        /*NOTE:Might actually work without the standard simplex*/
710                        dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
711                        for (int jj=1;jj<=this->numVars;jj++)
712                        {
713                                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
714                        }
715                        dd_WriteMatrix(stdout,intPointMatrix);
716                        intvec *iv_weight = new intvec(this->numVars);
717                        interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
718                        dd_FreeMatrix(intPointMatrix);
719                        dd_free_global_constants();
720                       
721                        /*Step 3
722                        turn the minimal basis into a reduced one
723                        */
724                        ring dstRing=rCopy0(srcRing);
725                        i=rBlocks(srcRing);
726                       
727                        dstRing->order=(int *)omAlloc((i+1)*sizeof(int));
728                        for(j=i;j>0;j--)
729                        {
730                                dstRing->order[j]=srcRing->order[j-1];
731                                dstRing->block0[j]=srcRing->block0[j-1];
732                                dstRing->block1[j]=srcRing->block1[j-1];
733                                if (srcRing->wvhdl[j-1] != NULL)
734                                {
735                                        dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
736                                }
737                        }
738                        dstRing->order[0]=ringorder_a;
739                        dstRing->order[1]=ringorder_dp;
740                        dstRing->order[2]=ringorder_C;                 
741                        dstRing->wvhdl[0] =( int *)omAlloc((iv_weight->length())*sizeof(int));
742                       
743                        for (int ii=0;ii<this->numVars;ii++)
744                        {                               
745                                dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
746                        }
747                        rComplete(dstRing);
748                        rChangeCurrRing(dstRing);
749#ifdef gfan_DEBUG
750                        rWrite(dstRing); cout << endl;
751#endif
752                        ideal dstRing_I;                       
753                        dstRing_I=idrCopyR(srcRing_HH,srcRing);                 
754                        //validOpts<1>=TRUE;
755#ifdef gfan_DEBUG
756                        idShow(dstRing_I);
757#endif                 
758                        BITSET save=test;
759                        test|=Sy_bit(OPT_REDSB);
760                        test|=Sy_bit(6);        //OPT_DEBUG                                     
761                        dstRing_I=kStd(idrCopyR(this->inputIdeal,this->rootRing),NULL,testHomog,NULL);                                 
762                        kInterRed(dstRing_I);
763                        idSkipZeroes(dstRing_I);
764                        test=save;
765                        /*End of step 3 - reduction*/
766                       
767                        f->setFlipGB(dstRing_I);//store the flipped GB
768#ifdef gfan_DEBUG
769                        cout << "Flipped GB is: " << endl;
770                        f->printFlipGB();
771#endif                 
772                }//void flip(ideal gb, facet *f)
773                               
774                /** \brief Compute the remainder of a polynomial by a given ideal
775                *
776                * Compute \f$ f^{\mathcal{G}} \f$
777                * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
778                * However, since we are only interested in the remainder, there is no need to
779                * compute the factors \f$ a_i \f$
780                */
781                //NOTE: Should be replaced by kNF or kNF2
782                poly restOfDiv(poly const &f, ideal const &I)
783                {
784                        cout << "Entering restOfDiv" << endl;
785                        poly p=f;
786                        pWrite(p);
787                        //poly r=kCreateZeroPoly(,currRing,currRing);   //The 0-polynomial, hopefully
788                        poly r=NULL;    //The zero polynomial
789                        int ii;
790                        bool divOccured;
791                       
792                        while (p!=NULL)
793                        {
794                                ii=1;
795                                divOccured=FALSE;
796                               
797                                while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
798                                {                                       
799                                        if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
800                                        {
801                                                //cout << "TICK 3" << endl;
802                                                poly step1,step2,step3;
803                                                //cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
804                                                step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
805                                                //cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
806                                                //cout << "TICK 3.1" << endl;
807                                                step2 = ppMult_qq(step1, I->m[ii-1]);
808                                                //cout << "TICK 3.2" << endl;
809                                                step3 = pSub(pCopy(p), step2);
810                                                //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
811                                                //cout << "TICK 4" << endl;
812                                                //pSetm(p);
813                                                pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
814                                                p=step3;
815                                                pWrite(p);                                             
816                                                divOccured=TRUE;
817                                        }
818                                        else
819                                        {
820                                                ii += 1;
821                                        }//if (pLmDivisibleBy(I->m[ii],p,currRing))
822                                }//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
823                                if (divOccured==FALSE)
824                                {
825                                        //cout << "TICK 5" << endl;
826                                        r=pAdd(pCopy(r),pHead(p));
827                                        pSetm(r);
828                                        pSort(r);
829                                        //cout << "r="; pWrite(r); cout << endl;
830                                        //cout << "TICK 6" << endl;
831                                        if (pLength(p)!=1)
832                                        {
833                                                p=pSub(pCopy(p),pHead(p));      //Here it may occur that p=0 instead of p=NULL
834                                        }
835                                        else
836                                        {
837                                                p=NULL; //Hack to correct this situation                                               
838                                        }
839                                        //cout << "TICK 7" << endl;
840                                        //cout << "p="; pWrite(p);
841                                }//if (divOccured==FALSE)
842                        }//while (p!=0)
843                        return r;
844                }//poly restOfDiv(poly const &f, ideal const &I)
845               
846                /** \brief Compute \f$ f-f^{\mathcal{G}} \f$
847                */
848                //NOTE: use kNF or kNF2 instead of restOfDivision
849                ideal ffG(ideal const &H, ideal const &G)
850                {
851                        cout << "Entering ffG" << endl;
852                        int size=IDELEMS(H);
853                        ideal res=idInit(size,1);
854                        poly temp1, temp2, temp3;       //polys to temporarily store values for pSub
855                        for (int ii=0;ii<size;ii++)
856                        {
857                                res->m[ii]=restOfDiv(H->m[ii],G);
858                                //res->m[ii]=kNF(H->m[ii],G);
859                                temp1=H->m[ii];
860                                temp2=res->m[ii];                               
861                                temp3=pSub(temp1, temp2);
862                                res->m[ii]=temp3;
863                                //res->m[ii]=pSub(temp1,temp2); //buggy
864                                //pSort(res->m[ii]);
865                                //pSetm(res->m[ii]);
866                                cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                         
867                        }                       
868                        return res;
869                }
870               
871                /** \brief Compute a Groebner Basis
872                *
873                * Computes the Groebner basis and stores the result in gcone::gcBasis
874                *\param ideal
875                *\return void
876                */
877                void getGB(ideal const &inputIdeal)
878                {
879                        ideal gb;
880                        gb=kStd(inputIdeal,NULL,testHomog,NULL);
881                        idSkipZeroes(gb);
882                        this->gcBasis=gb;       //write the GB into gcBasis
883                }//void getGB
884               
885                /** \brief The Generic Groebner Walk due to FJLT
886                * Needed for computing the search facet
887                */
888                ideal GenGrbWlk(ideal, ideal)
889                {
890                }//GGW
891
892
893                /** \brief Compute the dot product of two intvecs
894                *
895                */
896                int dotProduct(intvec const &iva, intvec const &ivb)                           
897                {
898                        //intvec iva=a;
899                        //intvec ivb=b;
900                        int res=0;
901                        for (int i=0;i<this->numVars;i++)
902                        {
903                                res = res+(iva[i]*ivb[i]);
904                        }
905                        return res;
906                }//int dotProduct
907
908                /** \brief Check whether two intvecs are parallel
909                *
910                * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
911                */
912                bool isParallel(intvec const &a, intvec const &b)
913                {                       
914                        int lhs,rhs;
915                        lhs=dotProduct(a,b)*dotProduct(a,b);
916                        rhs=dotProduct(a,a)*dotProduct(b,b);
917                        cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
918                        if (lhs==rhs)
919                        {
920                                return TRUE;
921                        }
922                        else
923                        {
924                                return FALSE;
925                        }
926                }//bool isParallel
927               
928                /** \brief Compute an interior point of a given cone
929                */
930                void interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
931                {
932                        dd_LPPtr lp,lpInt;
933                        dd_ErrorType err=dd_NoError;
934                        dd_LPSolverType solver=dd_DualSimplex;
935                        dd_LPSolutionPtr lpSol=NULL;
936                        dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
937                        dd_rowindex ddnewpos;
938                        dd_NumberType numb;     
939                        //M->representation=dd_Inequality;
940                        //M->objective-dd_LPMin;  //Not sure whether this is needed
941                        dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
942                        //cout << "TICK 1" << endl;
943                                               
944                        //dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
945                        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
946                        //cout << "Tick 2" << endl;
947                        //dd_WriteMatrix(stdout,M);
948                       
949                        lp=dd_Matrix2LP(M, &err);
950                        if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}                     
951                        if (lp==NULL){cout << "LP is NULL" << endl;}
952                        dd_WriteLP(stdout,lp);
953                        //cout << "Tick 3" << endl;
954                                               
955                        lpInt=dd_MakeLPforInteriorFinding(lp);
956                        if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
957                        dd_WriteLP(stdout,lpInt);
958                        //cout << "Tick 4" << endl;
959                       
960                        dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
961                        if (err!=dd_NoError)
962                        {
963                                cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
964                                dd_WriteErrorMessages(stdout, err);
965                        }
966                       
967                        //dd_LPSolve(lpInt,solver,&err);        //This will not result in a point from the relative interior
968                        if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
969                        //cout << "Tick 5" << endl;
970                                                                       
971                        //lpSol=dd_CopyLPSolution(lpInt);
972                        if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}                 
973                        //cout << "Tick 6" << endl;
974#ifdef gfan_DEBUG
975                        cout << "Interior point: ";
976#endif
977                        for (int ii=1; ii<(lpSol->d)-1;ii++)
978                        {
979#ifdef gfan_DEBUG
980                                dd_WriteNumber(stdout,lpSol->sol[ii]);
981#endif
982                                /* NOTE This works only as long as gmp returns fractions with the same denominator*/
983                                (iv)[ii-1]=(int)mpz_get_d(mpq_numref(lpSol->sol[ii]));  //looks evil, but does the trick
984                        }
985                        dd_FreeLPSolution(lpSol);
986                        dd_FreeLPData(lpInt);
987                        dd_FreeLPData(lp);
988                        set_free(ddlinset);
989                        set_free(ddredrows);                   
990                       
991                }//void interiorPoint(dd_MatrixPtr const &M)
992               
993                ring rCopyAndChangeWeight(ring const r, intvec const *ivw)                             
994                {
995                        ring res=rCopy0(r, FALSE, FALSE);
996                        int i=rBlocks(r);
997                        int j;
998
999                        res->order=(int *)omAlloc((i+1)*sizeof(int));
1000                        /*res->block0=(int *)omAlloc0((i+1)*sizeof(int));
1001                        res->block1=(int *)omAlloc0((i+1)*sizeof(int));
1002                        int ** wvhdl =(int **)omAlloc0((i+1)*sizeof(int**));*/
1003                        for(j=i;j>0;j--)
1004                        {
1005                                res->order[j]=r->order[j-1];
1006                                res->block0[j]=r->block0[j-1];
1007                                res->block1[j]=r->block1[j-1];
1008                                if (r->wvhdl[j-1] != NULL)
1009                                {
1010                                        res->wvhdl[j] = (int*) omMemDup(r->wvhdl[j-1]);
1011                                }
1012                        }
1013                        res->order[0]=ringorder_a;
1014                        res->order[1]=ringorder_dp;
1015                        res->order[2]=ringorder_C;
1016                        //res->wvhdl = wvhdl;
1017                       
1018                        res->wvhdl[0] =( int *)omAlloc((ivw->length())*sizeof(int));
1019                        for (int ii=0;ii<this->numVars;ii++)
1020                        {                               
1021                                res->wvhdl[0][ii]=(*ivw)[ii];                           
1022                        }                       
1023                       
1024                        rComplete(res);
1025                        return res;
1026                }//rCopyAndChange
1027               
1028                /**
1029                * Determines whether a given facet of a cone is the search facet of a neighbouring cone
1030                * This is done in the following way:
1031                * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1032                * that is first crossed during the generic walk.
1033                * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1034                * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1035                */
1036                bool isSearchFacet(gcone &tmpcone, facet &testfacet)
1037                {                               
1038                        ring actRing=currRing;
1039                        facet *facetPtr=(facet*)tmpcone.facetPtr;                       
1040                        facet *fMin=new facet(*facetPtr);       //Pointer to the "minimal" facet
1041                        //facet *fMin = new facet(tmpcone.facetPtr);
1042                        //fMin=tmpcone.facetPtr;                //Initialise to first facet of tmpcone
1043                        facet *fAct;
1044                        //fMin = testfacet;
1045                        //fMin->next=testfacet->next;
1046                        fAct = fMin;
1047                       
1048                        cout << endl<< fMin->next << endl;
1049                        cout << tmpcone.facetPtr->next << endl;
1050                        tmpcone.facetPtr->printNormal();
1051                        fAct=tmpcone.facetPtr->next;
1052                        fAct->printNormal();
1053                        /*fMin->printNormal();
1054                        fMin->next->printNormal();*/
1055                       
1056                        rChangeCurrRing(this->rootRing);        //because we compare the monomials in the rootring                     
1057                        poly p=NULL;
1058                        poly q=NULL;
1059                        intvec *p_weight = new intvec(this->numVars);                   
1060                        intvec *q_weight = new intvec(this->numVars);
1061                        intvec *sigma = new intvec(this->numVars);
1062                        sigma=tmpcone.getIntPoint();
1063                        cout << "pling" << endl;
1064                        int *u=(int *)omAlloc((this->numVars+1)*sizeof(int));
1065                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
1066                        int weight1,weight2;
1067                        while(fAct->next!=NULL) //ersetzen durch fAct
1068                        {
1069                                /* Get alpha_i and alpha_{i+1} */
1070                                p_weight=fMin->getFacetNormal();
1071                                q_weight=fMin->next->getFacetNormal();
1072                                p_weight->show(); 
1073                                q_weight->show();
1074                                /*Compute the dot product of sigma and alpha_{i,j}*/
1075                                weight1=dotProduct(sigma,p_weight);
1076                                weight2=dotProduct(sigma,q_weight);
1077                                cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1078                                /*Adjust alpha_i and alpha_i+1 accordingly*/
1079                                for(int ii=1;ii<=this->numVars;ii++)
1080                                {
1081                                        /*p_weight[ii]=weight1*(*p_weight)[ii];
1082                                        q_weight[ii]=weight2*(*q_weight)[ii];*/
1083                                        u[ii]=weight1*(*p_weight)[ii];
1084                                        v[ii]=weight2*(*q_weight)[ii];
1085                                }
1086                                cout << "PLONK" << endl;
1087                               
1088                                /*Now p_weight and q_weight need to be compared as exponent vectors*/
1089                                pSetExpV(p,u);
1090                                pSetExpV(q,v);
1091                                cout << "AARGH" << endl;
1092                                /*We want to check whether x^p < x^q
1093                                => want to check for return value 1 */
1094                                if (pLmCmp(p,q)==1)     //i.e. x^q is smaller
1095                                {
1096                                        fMin=fMin->next;
1097                                        fAct=fMin;
1098                                }
1099                                else
1100                                {
1101                                        fAct=fAct->next;
1102                                }
1103                                fAct=fAct->next;
1104                        }//while(tmpcone.facetPtr->next!=NULL)
1105                       
1106                        /*If testfacet was minimal then fMin should still point there */
1107                        //if (isParallel(fMin->getFacetNormal,testfacet.getFacetNormal))
1108                        if (fMin==tmpcone.facetPtr)
1109                        {
1110                                rChangeCurrRing(actRing);
1111                                return TRUE;
1112                        }
1113                        else 
1114                        {
1115                                rChangeCurrRing(actRing);
1116                                return FALSE;
1117                        }
1118                }//bool isSearchFacet
1119               
1120                void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1121                {
1122                        facet *fAct=new facet();
1123                        fAct = gcAct->facetPtr;                 
1124                       
1125                        while(fAct->next!=NULL)  //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
1126                        {
1127                                cout << "==========================================================================================="<< endl;
1128                                gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1129                                gcone *gcTmp = new gcone(*gcAct);
1130                                idShow(gcTmp->gcBasis);
1131                                gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1132#ifdef gfan_DEBUG
1133                                facet *f = new facet();
1134                                f=gcTmp->facetPtr;
1135                                while(f->next!=NULL)
1136                                {
1137                                        f->printNormal();
1138                                        f=f->next;                                     
1139                                }
1140#endif
1141                                gcTmp->showIntPoint();
1142                                /*recursive part goes gere*/
1143                                if (isSearchFacet(*gcTmp,(facet&)gcAct->facetPtr))
1144                                {
1145                                        gcAct->next=gcTmp;
1146                                        cout << "PING"<< endl;
1147                                        reverseSearch(gcTmp);
1148                                }
1149                                else
1150                                {
1151                                        delete gcTmp;
1152                                        /*NOTE remove fAct from linked list. It's no longer needed*/
1153                                }
1154                                /*recursion ends*/
1155                                fAct = fAct->next;             
1156                        }//while(fAct->next!=NULL)
1157                }//reverseSearch
1158};//class gcone
1159
1160ideal gfan(ideal inputIdeal)
1161{
1162        int numvar = pVariables; 
1163       
1164#ifdef gfan_DEBUG
1165        cout << "Now in subroutine gfan" << endl;
1166#endif
1167        ring inputRing=currRing;        // The ring the user entered
1168        ring rootRing;                  // The ring associated to the target ordering
1169        ideal res;     
1170        facet *fRoot;
1171       
1172        /*
1173        1. Select target order, say dp.
1174        2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
1175        3. getConeNormals
1176        */
1177       
1178        /* Construct a new ring which will serve as our root
1179        Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
1180        resolved 07.04.2009 MM
1181        */
1182        rootRing=rCopy0(currRing);
1183        rootRing->order[0]=ringorder_lp;
1184        //NOTE: Build ring accordiing to rCopyAndChangeWeight
1185        /*rootRing->order[0]=ringorder_a;
1186        rootRing->order[1]=ringorder_lp;
1187        rootRing->wvhdl[0] =( int *)omAlloc(numvar*sizeof(int));
1188        rootRing->wvhdl[0][1]=1;
1189        rootRing->wvhdl[0][2]=1;*/
1190        rComplete(rootRing);
1191        rChangeCurrRing(rootRing);
1192       
1193        /* Fetch the inputIdeal into our rootRing */
1194        map theMap=(map)idMaxIdeal(1);  //evil hack!
1195        theMap->preimage=NULL;  //neccessary?
1196        ideal rootIdeal;
1197        rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
1198#ifdef gfan_DEBUG
1199        cout << "Root ideal is " << endl;
1200        idShow(rootIdeal);
1201        cout << "The root ring is " << endl;
1202        rWrite(rootRing);
1203        cout << endl;
1204#endif 
1205       
1206        //gcone *gcRoot = new gcone();  //Instantiate the sink
1207        gcone *gcRoot = new gcone(rootRing,rootIdeal);
1208        gcone *gcAct;
1209        gcAct = gcRoot;
1210        gcAct->numVars=pVariables;
1211        gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
1212        idShow(gcAct->gcBasis);
1213        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
1214        //gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1215        /*Now it is time to compute the search facets, respectively start the reverse search.
1216        But since we are in the root all facets should be search facets. IS THIS TRUE?
1217        NOTE: Check for flippability is not very sophisticated
1218        */
1219        /*facet *fAct=new facet();
1220        fAct=gcAct->facetPtr;
1221        while(fAct->next!=NULL)
1222        {
1223                gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1224                gcone *gcTmp = new gcone(*gcAct);
1225                idShow(gcTmp->gcBasis);
1226                gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1227                gcTmp->getIntPoint();
1228                fAct = fAct->next;             
1229        }*/
1230        gcAct->reverseSearch(gcAct);
1231        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
1232        The return type will then be of type LIST_CMD
1233        Assume gfan has finished, thus we have enumerated all the cones
1234        Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
1235        Groebner Basis and merge this somehow with LIST_CMD
1236        => Count the cones!
1237        */
1238        rChangeCurrRing(rootRing);
1239        //res=gcAct->gcBasis;
1240        res=gcRoot->gcBasis;   
1241        return res;
1242        //return GBlist;
1243}
1244/*
1245Since gfan.cc is #included from extra.cc there must not be a int main(){} here
1246*/
1247#endif
Note: See TracBrowser for help on using the repository browser.