source: git/kernel/gfan.cc @ 57bfa2

spielwiese
Last change on this file since 57bfa2 was ecc6a3, checked in by Martin Monerjan, 15 years ago
contd work on ffG and restOfDiv. Tailred already included by Alg from CLO? git-svn-id: file:///usr/local/Singular/svn/trunk@11690 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.2 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009-04-14 11:05:52 $
5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.32 2009-04-14 11:05:52 monerjan Exp $
6$Id: gfan.cc,v 1.32 2009-04-14 11:05:52 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
25//Hacks for different working places
26#define ITWM
27
28#ifdef UNI
29#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
30#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
31#endif
32
33#ifdef HOME
34#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
35#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
36#endif
37
38#ifdef ITWM
39#include "/u/slg/monerjan/cddlib/include/setoper.h"
40#include "/u/slg/monerjan/cddlib/include/cdd.h"
41#endif
42
43#ifndef gfan_DEBUG
44#define gfan_DEBUG
45#endif
46
47//#include gcone.h
48
49/**
50*\brief Class facet
51*       Implements the facet structure as a linked list
52*
53*/
54class facet
55{
56        private:
57                /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
58                intvec *fNormal;
59               
60                /** \brief The Groebner basis on the other side of a shared facet
61                 *
62                 * In order not to have to compute the flipped GB twice we store the basis we already get
63                 * when identifying search facets. Thus in the next step of the reverse search we can
64                 * just copy the old cone and update the facet and the gcBasis.
65                 * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
66                 */
67                ideal flipGB;           //The Groebner Basis on the other side, computed via gcone::flip
68
69                       
70        public:
71                /** The default constructor. Do I need a constructor of type facet(intvec)? */
72                facet()                 
73                {
74                        // Pointer to next facet.  */
75                        /* Defaults to NULL. This way there is no need to check explicitly */
76                        this->next=NULL; 
77                }
78               
79                /** The default destructor */
80                ~facet(){;}
81               
82                /** Stores the facet normal \param intvec*/
83                void setFacetNormal(intvec *iv){
84                        this->fNormal = ivCopy(iv);
85                        //return;
86                }
87               
88                /** Hopefully returns the facet normal */
89                intvec *getFacetNormal()
90                {       
91                        //this->fNormal->show();        cout << endl;   
92                        return this->fNormal;
93                }
94               
95                /** Method to print the facet normal*/
96                void printNormal()
97                {
98                        fNormal->show();
99                }
100               
101                /** Store the flipped GB*/
102                void setFlipGB(ideal I)
103                {
104                        this->flipGB=I;
105                }
106               
107                /** Print the flipped GB*/
108                void printFlipGB()
109                {
110                        idShow(this->flipGB);
111                }
112               
113                bool isFlippable;       //flippable facet? Want to have cone->isflippable.facet[i]
114                bool isIncoming;        //Is the facet incoming or outgoing?
115                facet *next;            //Pointer to next facet
116};
117
118/**
119*\brief Implements the cone structure
120*
121* A cone is represented by a linked list of facet normals
122* @see facet
123*/
124/*class gcone
125finally this should become s.th. like gconelib.{h,cc} to provide an API
126*/
127class gcone
128{
129        private:
130                int numFacets;          //#of facets of the cone               
131
132        public:
133                /** \brief Default constructor.
134                *
135                * Initialises this->next=NULL and this->facetPtr=NULL
136                */
137                gcone()
138                {
139                        this->next=NULL;
140                        this->facetPtr=NULL;
141                }
142                ~gcone();               //destructor
143               
144                /** Pointer to the first facet */
145                facet *facetPtr;        //Will hold the adress of the first facet; set by gcone::getConeNormals
146                poly gcMarkedTerm;      //marked terms of the cone's Groebner basis
147                int numVars;            //#of variables in the ring
148               
149                /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
150                ideal gcBasis;          //GB of the cone, set by gcone::getGB();
151                gcone *next;            //Pointer to *previous* cone in search tree
152               
153                /** \brief Compute the normals of the cone
154                *
155                * This method computes a representation of the cone in terms of facet normals. It takes an ideal
156                * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
157                * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
158                * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
159                * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
160                * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
161                * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
162                */
163                void getConeNormals(const ideal &I)
164                {
165#ifdef gfan_DEBUG
166                        cout << "*** Computing Inequalities... ***" << endl;
167#endif         
168                        //All variables go here - except ineq matrix and *v, see below
169                        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
170                        int pCompCount;                 // # of terms in a poly
171                        poly aktpoly;
172                        int numvar = pVariables;        // # of variables in a polynomial (or ring?)
173                        int leadexp[numvar];            // dirty hack of exp.vects
174                        int aktexp[numvar];
175                        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
176                        dd_rowrange ddrows;
177                        dd_colrange ddcols;
178                        dd_rowset ddredrows;            // # of redundant rows in ddineq
179                        dd_rowset ddlinset;             // the opposite
180                        dd_rowindex ddnewpos;           // all to make dd_Canonicalize happy
181                        dd_NumberType ddnumb=dd_Real;   //Number type
182                        dd_ErrorType dderr=dd_NoError;  //
183                        // End of var declaration
184#ifdef gfan_DEBUG
185                        cout << "The Groebner basis has " << lengthGB << " elements" << endl;
186                        cout << "The current ring has " << numvar << " variables" << endl;
187#endif
188                        cols = numvar;
189               
190                        //Compute the # inequalities i.e. rows of the matrix
191                        rows=0; //Initialization
192                        for (int ii=0;ii<IDELEMS(I);ii++)
193                        {
194                                aktpoly=(poly)I->m[ii];
195                                rows=rows+pLength(aktpoly)-1;
196                        }
197#ifdef gfan_DEBUG
198                        cout << "rows=" << rows << endl;
199                        cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
200#endif
201                        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
202                        dd_set_global_constants();
203                        ddrows=rows;
204                        ddcols=cols;
205                        dd_MatrixPtr ddineq;            //Matrix to store the inequalities
206                        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
207               
208                        // We loop through each g\in GB and compute the resulting inequalities
209                        for (int i=0; i<IDELEMS(I); i++)
210                        {
211                                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
212                                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
213                                cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
214               
215                                int *v=(int *)omAlloc((numvar+1)*sizeof(int));
216                                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
217                               
218                                //Store leadexp for aktpoly
219                                for (int kk=0;kk<numvar;kk++)
220                                {
221                                        leadexp[kk]=v[kk+1];
222                                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
223                                        //but compute the diff below
224                                }
225               
226                               
227                                while (pNext(aktpoly)!=NULL) //move to next term until NULL
228                                {
229                                        aktpoly=pNext(aktpoly);
230                                        pSetm(aktpoly);         //doesn't seem to help anything
231                                        pGetExpV(aktpoly,v);
232                                        for (int kk=0;kk<numvar;kk++)
233                                        {
234                                                aktexp[kk]=v[kk+1];
235                                                //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
236                                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
237                                        }
238                                        aktmatrixrow=aktmatrixrow+1;
239                                }//while
240               
241                        } //for
242               
243                        //Maybe add another row to contain the constraints of the standard simplex?
244
245#ifdef gfan_DEBUG
246                        cout << "The inequality matrix is" << endl;
247                        dd_WriteMatrix(stdout, ddineq);
248#endif
249
250                        // The inequalities are now stored in ddineq
251                        // Next we check for superflous rows
252                        ddredrows = dd_RedundantRows(ddineq, &dderr);
253                        if (dderr!=dd_NoError)                  // did an error occur?
254                        {
255                                dd_WriteErrorMessages(stderr,dderr);    //if so tell us
256                        } else
257                        {
258                                cout << "Redundant rows: ";
259                                set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
260                        }//if dd_Error
261               
262                        //Remove reduntant rows here!
263                        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
264                        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
265                        ddcols = ddineq->colsize;
266#ifdef gfan_DEBUG
267                        cout << "Having removed redundancies, the normals now read:" << endl;
268                        dd_WriteMatrix(stdout,ddineq);
269                        cout << "Rows = " << ddrows << endl;
270                        cout << "Cols = " << ddcols << endl;
271#endif
272                       
273                        /*Write the normals into class facet*/
274#ifdef gfan_DEBUG
275                        cout << "Creating list of normals" << endl;
276#endif
277                        /*The pointer *fRoot should be the return value of this function*/
278                        facet *fRoot = new facet();     //instantiate new facet
279                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
280                        facet *fAct;                    //instantiate pointer to active facet
281                        fAct = fRoot;                   //Seems to do the trick. fRoot and fAct have to point to the same adress!
282                        std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
283                        for (int kk = 0; kk<ddrows; kk++)
284                        {
285                                intvec *load = new intvec(numvar);      //intvec to store a single facet normal that will then be stored via setFacetNormal
286                                for (int jj = 1; jj <ddcols; jj++)
287                                {
288                                        double *foo;
289                                        foo = (double*)ddineq->matrix[kk][jj];  //get entry from actual position
290#ifdef gfan_DEBUG
291                                        std::cout << "fAct is " << *foo << " at " << fAct << std::endl;
292#endif 
293                                        (*load)[jj-1] = (int)*foo;              //store typecasted entry at pos jj-1 of load
294                                        //check for flipability here
295                                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
296                                        {
297                                                //fAct->next = new facet();     //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
298                                        }
299                                }//for (int jj = 1; jj <ddcols; jj++)
300                                /*Quick'n'dirty hack for flippability*/ 
301                                bool isFlippable;                       
302                                for (int jj = 0; jj<this->numVars; jj++)
303                                {                                       
304                                        intvec *ivCanonical = new intvec(this->numVars);
305                                        (*ivCanonical)[jj]=1;                                   
306                                        if (dotProduct(load,ivCanonical)>=0)
307                                        {
308                                                isFlippable=FALSE;                                             
309                                        }
310                                        else
311                                        {
312                                                isFlippable=TRUE;
313                                        }                                       
314                                        delete ivCanonical;
315                                }//for (int jj = 0; jj<this->numVars; jj++)
316                                if (isFlippable==FALSE)
317                                {
318                                        cout << "Ignoring facet";
319                                        load->show();
320                                        //fAct->next=NULL;
321                                }
322                                else
323                                {       /*Now load should be full and we can call setFacetNormal*/
324                                        fAct->setFacetNormal(load);
325                                        fAct->next = new facet();
326                                        //fAct->printNormal();
327                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
328                                }//if (isFlippable==FALSE)
329                                delete load;
330                        }//for (int kk = 0; kk<ddrows; kk++)
331                        /*
332                        Now we should have a linked list containing the facet normals of those facets that are
333                        -irredundant
334                        -flipable
335                        Adressing is done via *facetPtr
336                        */
337                       
338                        //Clean up but don't delete the return value! (Whatever it will turn out to be)
339                        dd_FreeMatrix(ddineq);
340                        set_free(ddredrows);
341                        free(ddnewpos);
342                        set_free(ddlinset);
343                        dd_free_global_constants();
344
345                }//method getConeNormals(ideal I)       
346               
347               
348                /** \brief Compute the Groebner Basis on the other side of a shared facet
349                *
350                * Implements algorithm 4.3.2 from Anders' thesis.
351                * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
352                * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
353                * is parallel to \f$ leadexp(g) \f$
354                * Parallelity is checked using basic linear algebra. See gcone::isParallel.
355                * Other possibilities includes  computing the rank of the matrix consisting of the vectors in question and
356                * computing an interior point of the facet and taking all terms having the same weight with respect
357                * to this interior point.
358                *\param ideal, facet
359                * Input: a marked,reduced Groebner basis and a facet
360                */
361                void flip(ideal gb, facet *f)           //Compute "the other side"
362                {                       
363                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
364                        fNormal = f->getFacetNormal();  //read this->fNormal;
365#ifdef gfan_DEBUG
366                        cout << "===" << endl;
367                        cout << "running gcone::flip" << endl;
368                        cout << "fNormal=";
369                        fNormal->show();
370                        cout << endl;
371#endif                         
372                        /*1st step: Compute the initial ideal*/
373                        poly initialFormElement[IDELEMS(gb)];   //array of #polys in GB to store initial form
374                        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
375                        poly aktpoly;
376                        intvec *check = new intvec(this->numVars);      //array to store the difference of LE and v
377                       
378                        for (int ii=0;ii<IDELEMS(gb);ii++)
379                        {
380                                aktpoly = (poly)gb->m[ii];                                                             
381                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
382                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
383                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
384                                initialFormElement[ii]=pHead(aktpoly);
385                               
386                                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
387                                {
388                                        aktpoly=pNext(aktpoly); //next term
389                                        pSetm(aktpoly);
390                                        pGetExpV(aktpoly,v);           
391                                        /* Convert (int)v into (intvec)check */                 
392                                        for (int jj=0;jj<this->numVars;jj++)
393                                        {
394                                                //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
395                                                //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
396                                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
397                                        }
398#ifdef gfan_DEBUG
399                                        cout << "check=";                       
400                                        check->show();
401                                        cout << endl;
402#endif
403                                        if (isParallel(check,fNormal)) //pass *check when
404                                        {
405                                                cout << "Parallel vector found, adding to initialFormElement" << endl;                         
406                                                initialFormElement[ii] = pAdd(initialFormElement[ii],(poly)pHead(aktpoly));
407                                        }                                               
408                                }//while
409                                cout << "Initial Form=";                               
410                                pWrite(initialFormElement[ii]);
411                                cout << "---" << endl;
412                                /*Now initialFormElement must be added to (ideal)initialForm */
413                                initialForm->m[ii]=initialFormElement[ii];
414                        }//for
415                        f->setFlipGB(initialForm);                     
416#ifdef gfan_DEBUG
417                        cout << "Initial ideal is: " << endl;
418                        idShow(initialForm);
419                        f->printFlipGB();
420                        cout << "===" << endl;
421#endif
422                        delete check;
423                       
424                        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
425                        /*Substep 2.1
426                        compute $G_{-\alpha}(in_v(I))
427                        see journal p. 66
428                        */
429                        ring srcRing=currRing;
430                        ring dstRing=rCopy0(srcRing);
431                        dstRing->order[0]=ringorder_a;
432                        //tmpring->order[1]=ringorder_dp;
433                        //tmpring->order[2]=ringorder_C;
434                        dstRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
435                       
436                        for (int ii=0;ii<this->numVars;ii++)
437                        {                               
438                                dstRing->wvhdl[0][ii]=-(*fNormal)[ii];  //What exactly am I doing here?
439                                //cout << tmpring->wvhdl[0][ii] << endl;
440                        }
441                        rComplete(dstRing);
442                        rChangeCurrRing(dstRing);
443                        //map theMap=(map)idMaxIdeal(1);
444                        rWrite(currRing); cout << endl;
445                        ideal ina;
446                        //ina=fast_map(initialForm,srcRing,(ideal)theMap,dstRing);                     
447                        ina=idrCopyR(initialForm,srcRing);
448                        cout << "ina=";
449                        idShow(ina); cout << endl;
450                        ideal H;
451                        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
452                        idSkipZeroes(H);
453                        cout << "H="; idShow(H); cout << endl;
454                        /*Substep 2.2
455                        do the lifting
456                        */
457                        rChangeCurrRing(srcRing);
458                        ideal srcRing_H;
459                        ideal srcRing_HH;
460                        //map theMap = (map)idMaxIdeal(1);
461                        //srcRing_H=fast_map(H,dstRing,(ideal)theMap,srcRing);
462                        srcRing_H=idrCopyR(H,dstRing);
463                        idShow(srcRing_H);
464                        srcRing_HH=ffG(srcRing_H,this->gcBasis);                       
465                        idShow(srcRing_HH);
466                        /*Substep 2.3
467                        turn the minimal basis into a reduced one
468                        */
469                       
470                }//void flip(ideal gb, facet *f)
471                               
472                /** \brief Compute the remainder of a polynomial by a given ideal
473                *
474                * Compute \f$ f^{\mathcal{G}} \f$
475                * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
476                * However, since we are only interested in the remainder, there is no need to
477                * compute the factors \f$ a_i \f$
478                */
479                poly restOfDiv(poly const &f, ideal const &I)
480                {
481                        cout << "Entering restOfDiv" << endl;
482                        poly p=f;
483                        pWrite(p);
484                        //poly r=kCreateZeroPoly(,currRing,currRing);   //The 0-polynomial, hopefully
485                        poly r=NULL;    //The zero polynomial
486                        int ii;
487                        bool divOccured;
488                       
489                        while (p!=NULL)
490                        {
491                                ii=1;
492                                divOccured=FALSE;
493                               
494                                while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
495                                {                                       
496                                        if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
497                                        {
498                                                //cout << "TICK 3" << endl;
499                                                poly step1,step2,step3;
500                                                cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
501                                                step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
502                                                cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
503                                                //cout << "TICK 3.1" << endl;
504                                                step2 = ppMult_qq(step1, I->m[ii-1]);
505                                                //cout << "TICK 3.2" << endl;
506                                                step3 = pSub(p, step2);
507                                                //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
508                                                //cout << "TICK 4" << endl;
509                                                pSetm(p);
510                                                p=step3;
511                                                pWrite(p);                                             
512                                                divOccured=TRUE;
513                                        }
514                                        else
515                                        {
516                                                ii += 1;
517                                        }//if (pLmDivisibleBy(I->m[ii],p,currRing))
518                                }//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
519                                if (divOccured==FALSE)
520                                {
521                                        //cout << "TICK 5" << endl;
522                                        r=pAdd(r,pHead(p));
523                                        pSetm(r);
524                                        cout << "r="; pWrite(r); cout << endl;
525                                        //cout << "TICK 6" << endl;
526                                        if (pLength(p)!=1)
527                                        {
528                                                p=pSub(p,pHead(p));     //Here it may occur that p=0 instead of p=NULL
529                                        }
530                                        else
531                                        {
532                                                p=NULL; //Hack to correct this situation                                               
533                                        }
534                                        //cout << "TICK 7" << endl;
535                                        cout << "p="; pWrite(p);
536                                }//if (divOccured==FALSE)
537                        }//while (p!=0)
538                        return r;
539                }//poly restOfDiv(poly const &f, ideal const &I)
540               
541                /** \brief Compute \f$ f-f^{\mathcal{G}} \f$
542                */
543                ideal ffG(ideal const &H, ideal const &G)
544                {
545                        cout << "Entering ffG" << endl;
546                        int size=IDELEMS(H);
547                        ideal res=idInit(size,1);
548                        for (int ii=0;ii<size;ii++)
549                        {
550                                res->m[ii]=restOfDiv(H->m[ii],G);
551                                //res->m[ii]=pSub(H->m[ii],res->m[ii]); //buggy
552                                pSetm(res->m[ii]);
553                                cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                         
554                        }
555                        return res;
556                }
557               
558                /** \brief Compute a Groebner Basis
559                *
560                * Computes the Groebner basis and stores the result in gcone::gcBasis
561                *\param ideal
562                *\return void
563                */
564                void getGB(ideal const &inputIdeal)
565                {
566                        ideal gb;
567                        gb=kStd(inputIdeal,NULL,testHomog,NULL);
568                        idSkipZeroes(gb);
569                        this->gcBasis=gb;       //write the GB into gcBasis
570                }//void getGB
571               
572                ideal GenGrbWlk(ideal, ideal);  //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
573
574
575                /** \brief Compute the dot product of two intvecs
576                *
577                */
578                int dotProduct(intvec const &iva, intvec const &ivb)                           
579                {
580                        //intvec iva=a;
581                        //intvec ivb=b;
582                        int res=0;
583                        for (int i=0;i<this->numVars;i++)
584                        {
585                                res = res+(iva[i]*ivb[i]);
586                        }
587                        return res;
588                }//int dotProduct
589
590                /** \brief Check whether two intvecs are parallel
591                *
592                * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
593                */
594                bool isParallel(intvec const &a, intvec const &b)
595                {                       
596                        int lhs,rhs;
597                        lhs=dotProduct(a,b)*dotProduct(a,b);
598                        rhs=dotProduct(a,a)*dotProduct(b,b);
599                        cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
600                        if (lhs==rhs)
601                        {
602                                return TRUE;
603                        }
604                        else
605                        {
606                                return FALSE;
607                        }
608                }//bool isParallel
609};//class gcone
610
611ideal gfan(ideal inputIdeal)
612{
613        int numvar = pVariables; 
614       
615        #ifdef gfan_DEBUG
616        cout << "Now in subroutine gfan" << endl;
617        #endif
618        ring inputRing=currRing;        // The ring the user entered
619        ring rootRing;                  // The ring associated to the target ordering
620        ideal res;
621        //matrix ineq;                  //Matrix containing the boundary inequalities
622        facet *fRoot;
623       
624        /*
625        1. Select target order, say dp.
626        2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
627        3. getConeNormals
628        */
629       
630        /* Construct a new ring which will serve as our root
631        Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
632        resolved 07.04.2009 MM
633        */
634        rootRing=rCopy0(currRing);
635        rootRing->order[0]=ringorder_lp;
636        rComplete(rootRing);
637        rChangeCurrRing(rootRing);
638       
639        /* Fetch the inputIdeal into our rootRing */
640        map theMap=(map)idMaxIdeal(1);  //evil hack!
641        theMap->preimage=NULL;  //neccessary?
642        ideal rootIdeal;
643        rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
644#ifdef gfan_DEBUG
645        cout << "Root ideal is " << endl;
646        idShow(rootIdeal);
647        cout << "The root ring is " << endl;
648        rWrite(rootRing);
649        cout << endl;
650#endif 
651       
652        gcone *gcRoot = new gcone();    //Instantiate the sink
653        gcone *gcAct;
654        gcAct = gcRoot;
655        gcAct->numVars=pVariables;
656        gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
657        idShow(gcAct->gcBasis);
658        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
659        gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
660        /*Now it is time to compute the search facets, respectively start the reverse search.
661        But since we are in the root all facets should be search facets. IS THIS TRUE?
662        MIND: AS OF NOW, THE LIST OF FACETS IS NOT PURGED OF NON-FLIPPAPLE FACETS
663        */
664       
665        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
666        The return type will then be of type LIST_CMD
667        Assume gfan has finished, thus we have enumerated all the cones
668        Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
669        Groebner Basis and merge this somehow with LIST_CMD
670        => Count the cones!
671        */
672        rChangeCurrRing(inputRing);
673        res=gcAct->gcBasis;
674        //cout << fRoot << endl;
675        return res;
676        //return GBlist;
677}
678/*
679Since gfan.cc is #included from extra.cc there must not be a int main(){} here
680*/
681#endif
Note: See TracBrowser for help on using the repository browser.