Context Navigation

source: git/kernel/gfan.cc @ 459b80

spielwiese

Last change on this file since 459b80 was 459b80, checked in by Martin Monerjan, 14 years ago
Started work on markings git-svn-id: file:///usr/local/Singular/svn/trunk@11724 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 23.7 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009-04-16 09:59:16 $
5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.34 2009-04-16 09:59:16 monerjan Exp $
6	$Id: gfan.cc,v 1.34 2009-04-16 09:59:16 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	*/
54	class 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
125	finally this should become s.th. like gconelib.{h,cc} to provide an API
126	*/
127	class 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(ideal const &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(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
407	}
408	}//while
409	#ifdef gfan_DEBUG
410	cout << "Initial Form=";
411	pWrite(initialFormElement[ii]);
412	cout << "---" << endl;
413	#endif
414	/Now initialFormElement must be added to (ideal)initialForm /
415	initialForm->m[ii]=initialFormElement[ii];
416	}//for
417	f->setFlipGB(initialForm); //FIXME PROBABLY WRONG TO STORE HERE SINCE INA!=flibGB
418	#ifdef gfan_DEBUG
419	cout << "Initial ideal is: " << endl;
420	idShow(initialForm);
421	f->printFlipGB();
422	cout << "===" << endl;
423	#endif
424	delete check;
425
426	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
427	/*Substep 2.1
428	compute $G_{-\alpha}(in_v(I))
429	see journal p. 66
430	*/
431	ring srcRing=currRing;
432	ring dstRing=rCopy0(srcRing);
433	dstRing->order[0]=ringorder_a;
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
444	rWrite(currRing); cout << endl;
445	ideal ina;
446	ina=idrCopyR(initialForm,srcRing);
447	#ifdef gfan_DEBUG
448	cout << "ina=";
449	idShow(ina); cout << endl;
450	#endif
451	ideal H;
452	H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
453	idSkipZeroes(H);
454	#ifdef gfan_DEBUG
455	cout << "H="; idShow(H); cout << endl;
456	#endif
457	/*Substep 2.2
458	do the lifting and mark according to H
459	*/
460	rChangeCurrRing(srcRing);
461	ideal srcRing_H;
462	ideal srcRing_HH;
463	srcRing_H=idrCopyR(H,dstRing);
464	idShow(srcRing_H);
465	srcRing_HH=ffG(srcRing_H,this->gcBasis);
466	idShow(srcRing_HH);
467	/*Substep 2.2.1
468	Mark according to G_-\alpha
469	Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
470	we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
471	represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
472	Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
473	compute the difference accordingly
474	*/
475	dd_set_global_constants;
476	bool markingsAreCorrect=FALSE;
477	dd_MatrixPtr intPointMatrix;
478	int iPMatrixRows=0;
479	dd_rowrange aktrow=0;
480	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
481	{
482	poly aktpoly=(poly)srcRing_HH->m[ii];
483	iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
484	}
485	intPointMatrix = dd_CreateMatrix(iPMatrixRows,this->numVars+1);
486
487	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
488	{
489	poly aktpoly=srcRing_HH->m[ii];
490	/Comparison of leading monomials is done via exponent vectors/
491	for (int jj=0;jj<IDELEMS(H);jj++)
492	{
493	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
494	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
495	pGetExpV(aktpoly,src_ExpV);
496	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
497	cout << *src_ExpV << endl;
498	cout << *dst_ExpV << endl;
499	if (src_ExpV == dst_ExpV)
500	{
501	markingsAreCorrect=TRUE; //everything is fine
502	cout << "correct markings" << endl;
503	}//if (pHead(aktpoly)==pHead(H->m[jj])
504	delete src_ExpV;
505	delete dst_ExpV;
506	}//for (int jj=0;jj<IDELEMS(H);jj++)
507
508	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
509	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
510	if (markingsAreCorrect==TRUE)
511	{
512	pGetExpV(aktpoly,leadExpV);
513	}
514	else
515	{
516	pGetExpV(pCopy(pHead(H->m[ii])),leadExpV); //We use H->m[ii] as leading monomial
517	}
518	/compute differences of the expvects/
519	while (pNext(aktpoly)!=NULL)
520	{
521	aktpoly=pNext(aktpoly);
522	pGetExpV(aktpoly,v);
523	for (int jj=0;jj<this->numVars;jj++)
524	{
525	/Store into ddMatrix/
526	/FIXME Wrong values/
527	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],v[jj+1]-leadExpV[jj+1]);
528	}
529	aktrow +=1;
530	}
531	delete v;
532	delete leadExpV;
533	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
534	dd_WriteMatrix(stdout,intPointMatrix);
535	dd_FreeMatrix(intPointMatrix);
536
537	/*Step 3
538	turn the minimal basis into a reduced one
539	*/
540
541	}//void flip(ideal gb, facet *f)
542
543	/** \brief Compute the remainder of a polynomial by a given ideal
544	*
545	* Compute \f$ f^{\mathcal{G}} \f$
546	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
547	* However, since we are only interested in the remainder, there is no need to
548	* compute the factors \f$ a_i \f$
549	*/
550	poly restOfDiv(poly const &f, ideal const &I)
551	{
552	cout << "Entering restOfDiv" << endl;
553	poly p=f;
554	pWrite(p);
555	//poly r=kCreateZeroPoly(,currRing,currRing); //The 0-polynomial, hopefully
556	poly r=NULL; //The zero polynomial
557	int ii;
558	bool divOccured;
559
560	while (p!=NULL)
561	{
562	ii=1;
563	divOccured=FALSE;
564
565	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
566	{
567	if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
568	{
569	//cout << "TICK 3" << endl;
570	poly step1,step2,step3;
571	cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
572	step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
573	cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
574	//cout << "TICK 3.1" << endl;
575	step2 = ppMult_qq(step1, I->m[ii-1]);
576	//cout << "TICK 3.2" << endl;
577	step3 = pSub(pCopy(p), step2);
578	//p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
579	//cout << "TICK 4" << endl;
580	//pSetm(p);
581	pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
582	p=step3;
583	pWrite(p);
584	divOccured=TRUE;
585	}
586	else
587	{
588	ii += 1;
589	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
590	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
591	if (divOccured==FALSE)
592	{
593	//cout << "TICK 5" << endl;
594	r=pAdd(pCopy(r),pHead(p));
595	pSetm(r);
596	pSort(r);
597	cout << "r="; pWrite(r); cout << endl;
598	//cout << "TICK 6" << endl;
599	if (pLength(p)!=1)
600	{
601	p=pSub(pCopy(p),pHead(p)); //Here it may occur that p=0 instead of p=NULL
602	}
603	else
604	{
605	p=NULL; //Hack to correct this situation
606	}
607	//cout << "TICK 7" << endl;
608	cout << "p="; pWrite(p);
609	}//if (divOccured==FALSE)
610	}//while (p!=0)
611	return r;
612	}//poly restOfDiv(poly const &f, ideal const &I)
613
614	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
615	*/
616	ideal ffG(ideal const &H, ideal const &G)
617	{
618	cout << "Entering ffG" << endl;
619	int size=IDELEMS(H);
620	ideal res=idInit(size,1);
621	poly temp1, temp2, temp3; //polys to temporarily store values for pSub
622	for (int ii=0;ii<size;ii++)
623	{
624	res->m[ii]=restOfDiv(H->m[ii],G);
625	temp1=H->m[ii];
626	temp2=res->m[ii];
627	temp3=pSub(temp1, temp2);
628	res->m[ii]=temp3;
629	//res->m[ii]=pSub(temp1,temp2); //buggy
630	//pSort(res->m[ii]);
631	//pSetm(res->m[ii]);
632	cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
633	}
634	return res;
635	}
636
637	/** \brief Compute a Groebner Basis
638	*
639	* Computes the Groebner basis and stores the result in gcone::gcBasis
640	*\param ideal
641	*\return void
642	*/
643	void getGB(ideal const &inputIdeal)
644	{
645	ideal gb;
646	gb=kStd(inputIdeal,NULL,testHomog,NULL);
647	idSkipZeroes(gb);
648	this->gcBasis=gb; //write the GB into gcBasis
649	}//void getGB
650
651	ideal GenGrbWlk(ideal, ideal); //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
652
653
654	/** \brief Compute the dot product of two intvecs
655	*
656	*/
657	int dotProduct(intvec const &iva, intvec const &ivb)
658	{
659	//intvec iva=a;
660	//intvec ivb=b;
661	int res=0;
662	for (int i=0;i<this->numVars;i++)
663	{
664	res = res+(iva[i]*ivb[i]);
665	}
666	return res;
667	}//int dotProduct
668
669	/** \brief Check whether two intvecs are parallel
670	*
671	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
672	*/
673	bool isParallel(intvec const &a, intvec const &b)
674	{
675	int lhs,rhs;
676	lhs=dotProduct(a,b)*dotProduct(a,b);
677	rhs=dotProduct(a,a)*dotProduct(b,b);
678	cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
679	if (lhs==rhs)
680	{
681	return TRUE;
682	}
683	else
684	{
685	return FALSE;
686	}
687	}//bool isParallel
688	};//class gcone
689
690	ideal gfan(ideal inputIdeal)
691	{
692	int numvar = pVariables;
693
694	#ifdef gfan_DEBUG
695	cout << "Now in subroutine gfan" << endl;
696	#endif
697	ring inputRing=currRing; // The ring the user entered
698	ring rootRing; // The ring associated to the target ordering
699	ideal res;
700	//matrix ineq; //Matrix containing the boundary inequalities
701	facet *fRoot;
702
703	/*
704	1. Select target order, say dp.
705	2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
706	3. getConeNormals
707	*/
708
709	/* Construct a new ring which will serve as our root
710	Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
711	resolved 07.04.2009 MM
712	*/
713	rootRing=rCopy0(currRing);
714	rootRing->order[0]=ringorder_lp;
715	rComplete(rootRing);
716	rChangeCurrRing(rootRing);
717
718	/* Fetch the inputIdeal into our rootRing */
719	map theMap=(map)idMaxIdeal(1); //evil hack!
720	theMap->preimage=NULL; //neccessary?
721	ideal rootIdeal;
722	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
723	#ifdef gfan_DEBUG
724	cout << "Root ideal is " << endl;
725	idShow(rootIdeal);
726	cout << "The root ring is " << endl;
727	rWrite(rootRing);
728	cout << endl;
729	#endif
730
731	gcone *gcRoot = new gcone(); //Instantiate the sink
732	gcone *gcAct;
733	gcAct = gcRoot;
734	gcAct->numVars=pVariables;
735	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
736	idShow(gcAct->gcBasis);
737	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
738	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
739	/*Now it is time to compute the search facets, respectively start the reverse search.
740	But since we are in the root all facets should be search facets. IS THIS TRUE?
741	MIND: AS OF NOW, THE LIST OF FACETS IS NOT PURGED OF NON-FLIPPAPLE FACETS
742	*/
743
744	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
745	The return type will then be of type LIST_CMD
746	Assume gfan has finished, thus we have enumerated all the cones
747	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
748	Groebner Basis and merge this somehow with LIST_CMD
749	=> Count the cones!
750	*/
751	rChangeCurrRing(inputRing);
752	res=gcAct->gcBasis;
753	//cout << fRoot << endl;
754	return res;
755	//return GBlist;
756	}
757	/*
758	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
759	*/
760	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats: