Context Navigation

source: git/kernel/gfan.cc @ f26ea0

spielwiese

Last change on this file since f26ea0 was f26ea0, checked in by Martin Monerjan, 15 years ago
Maintenance and bugfixes method gcone::showFacets() now takes the codim as an optional parameter New method gcone::getCodim2Normals Made UCN global and modified constructor of gcone accordingly git-svn-id: file:///usr/local/Singular/svn/trunk@11888 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 51.0 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009-06-09 15:23:08 $
5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.58 2009-06-09 15:23:08 monerjan Exp $
6	$Id: gfan.cc,v 1.58 2009-06-09 15:23:08 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>
24	#include <bitset>
25	#include <fstream> //read-write cones to files
26	#include <gmp.h>
27	//#include <gmpxx.h>
28
29	/DO NOT REMOVE THIS/
30	#ifndef GMPRATIONAL
31	#define GMPRATIONAL
32	#endif
33
34	//Hacks for different working places
35	#define ITWM
36
37	#ifdef UNI
38	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
39	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
40	#endif
41
42	#ifdef HOME
43	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
44	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
45	#endif
46
47	#ifdef ITWM
48	#include "/u/slg/monerjan/cddlib/include/setoper.h"
49	#include "/u/slg/monerjan/cddlib/include/cdd.h"
50	#include "/u/slg/monerjan/cddlib/include/cddmp.h"
51	#endif
52
53	#ifndef gfan_DEBUG
54	#define gfan_DEBUG
55	#endif
56
57	//#include gcone.h
58	using namespace std;
59
60	/**
61	*\brief Class facet
62	* Implements the facet structure as a linked list
63	*
64	*/
65	class facet
66	{
67	private:
68	/** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
69	intvec *fNormal;
70
71	/** \brief An interior point of the facet*/
72	intvec *interiorPoint;
73
74	/** \brief Universal Cone Number
75	* The number of the cone the facet belongs to
76	*/
77	int UCN;
78
79	/** \brief The codim of the facet
80	*/
81	int codim;
82
83	/** \brief The Groebner basis on the other side of a shared facet
84	*
85	* In order not to have to compute the flipped GB twice we store the basis we already get
86	* when identifying search facets. Thus in the next step of the reverse search we can
87	* just copy the old cone and update the facet and the gcBasis.
88	* facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
89	*/
90	ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip
91
92	public:
93	//bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i]
94	bool isIncoming; //Is the facet incoming or outgoing?
95	facet *next; //Pointer to next facet
96	facet *codim2Ptr; //Pointer to (codim-2)-facet. Bit of recursion here ;-)
97	//intvec codim2Normals =(intvec)omAlloc0(sizeof(intvec*)); //Integer matrix containing the (codim-2)-facets
98
99	/** The default constructor. Do I need a constructor of type facet(intvec)? */
100	facet()
101	{
102	// Pointer to next facet. */
103	/* Defaults to NULL. This way there is no need to check explicitly */
104	this->next=NULL;
105	this->UCN=NULL;
106	this->codim2Ptr=NULL;
107	this->codim=1; //default for (codim-1)-facets
108	}
109
110	/** \brief Constructor for facets of codim >= 2
111	*/
112	facet(int const &n)
113	{
114	this->next=NULL;
115	this->UCN=NULL;
116	this->codim2Ptr=NULL;
117	if(n>1)
118	{
119	this->codim=n;
120	}//NOTE Handle exception here!
121	}
122
123	/** The default destructor */
124	~facet(){;}
125
126	/** Stores the facet normal \param intvec*/
127	void setFacetNormal(intvec *iv)
128	{
129	this->fNormal = ivCopy(iv);
130	}
131
132	/** Hopefully returns the facet normal */
133	intvec *getFacetNormal()
134	{
135	return this->fNormal;
136	}
137
138	/** Method to print the facet normal*/
139	void printNormal()
140	{
141	fNormal->show();
142	}
143
144	/** Store the flipped GB*/
145	void setFlipGB(ideal I)
146	{
147	this->flipGB=I;
148	}
149
150	/** Return the flipped GB*/
151	ideal getFlipGB()
152	{
153	return this->flipGB;
154	}
155
156	/** Print the flipped GB*/
157	void printFlipGB()
158	{
159	idShow(this->flipGB);
160	}
161
162	void setUCN(int n)
163	{
164	this->UCN=n;
165	}
166
167	int getUCN()
168	{
169	return this->UCN;
170	}
171
172	void setInteriorPoint(intvec *iv)
173	{
174	this->interiorPoint = ivCopy(iv);
175	}
176
177	intvec *getInteriorPoint()
178	{
179	return this->interiorPoint;
180	}
181
182	/*bool isFlippable(intvec &load)
183	{
184	bool res=TRUE;
185	int jj;
186	for (int jj = 0; jj<load.length(); jj++)
187	{
188	intvec *ivCanonical = new intvec(load.length());
189	(*ivCanonical)[jj]=1;
190	if (ivMult(&load,ivCanonical)<0)
191	{
192	res=FALSE;
193	break;
194	}
195	}
196	return res;
197
198	/*while (dotProduct(load,ivCanonical)>=0)
199	{
200	if (jj!=this->numVars)
201	{
202	intvec *ivCanonical = new intvec(this->numVars);
203	(*ivCanonical)[jj]=1;
204	res=TRUE;
205	jj += 1;
206	}
207	}
208	if (jj==this->numVars)
209	{
210	delete ivCanonical;
211	return FALSE;
212	}
213	else
214	{
215	delete ivCanonical;
216	return TRUE;
217	}*/
218	//}//bool isFlippable(facet &f)
219
220
221	friend class gcone; //Bad style
222	};
223
224	/**
225	*\brief Implements the cone structure
226	*
227	* A cone is represented by a linked list of facet normals
228	* @see facet
229	*/
230	/*class gcone
231	finally this should become s.th. like gconelib.{h,cc} to provide an API
232	*/
233	class gcone
234	{
235	private:
236	ring rootRing; //good to know this -> generic walk
237	ideal inputIdeal; //the original
238	ring baseRing; //the basering of the cone
239	/* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
240	intvec *ivIntPt; //an interior point of the cone
241	static int UCN; //unique number of the cone
242
243	public:
244	/** \brief Pointer to the first facet */
245	facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals
246
247	/** # of variables in the ring */
248	int numVars; //#of variables in the ring
249
250	/** # of facets of the cone
251	* This value is set by gcone::getConeNormals
252	*/
253	int numFacets; //#of facets of the cone
254
255	/** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
256	ideal gcBasis; //GB of the cone, set by gcone::getGB();
257	gcone next; //Pointer to previous* cone in search tree
258
259	/** \brief Default constructor.
260	*
261	* Initialises this->next=NULL and this->facetPtr=NULL
262	*/
263	gcone()
264	{
265	this->next=NULL;
266	this->facetPtr=NULL; //maybe this->facetPtr = new facet();
267	this->baseRing=currRing;
268	this->UCN++;
269	this->numFacets=0;
270	}
271
272	/** \brief Constructor with ring and ideal
273	*
274	* This constructor takes the root ring and the root ideal as parameters and stores
275	* them in the private members gcone::rootRing and gcone::inputIdeal
276	* Since knowledge of the root ring is only needed when using reverse search,
277	* this constructor is not needed when using the "second" method
278	*/
279	gcone(ring r, ideal I)
280	{
281	this->next=NULL;
282	this->facetPtr=NULL;
283	this->rootRing=r;
284	this->inputIdeal=I;
285	this->baseRing=currRing;
286	this->UCN++;
287	this->numFacets=0;
288	}
289
290	/** \brief Copy constructor
291	*
292	* Copies a cone, sets this->gcBasis to the flipped GB and reverses the
293	* direction of the according facet normal
294	* Call this only after a successfull call to gcone::flip which sets facet::flipGB
295	*/
296	//NOTE Prolly need to specify the facet to flip over
297	gcone(const gcone& gc, const facet &f)
298	{
299	this->next=NULL;
300	this->numVars=gc.numVars;
301	this->UCN=(gc.UCN)+1; //add 1 to the UCN of previous cone. This is NOT UNIQUE!
302	facet *fAct= new facet();
303	this->facetPtr=fAct;
304
305	intvec *ivtmp = new intvec(this->numVars);
306	//NOTE ivtmp = f->getFacetNormal();
307	ivtmp = gc.facetPtr->getFacetNormal();
308
309	ideal gb;
310	//NOTE gb=f->getFlipGB();
311	gb=gc.facetPtr->getFlipGB();
312	this->gcBasis=gb; //this cone's GB is the flipped GB
313
314	/Reverse direction of the facet normal to make it an inner normal/
315	for (int ii=0; ii<this->numVars;ii++)
316	{
317	(ivtmp)[ii]=-(ivtmp)[ii];
318	}
319
320	fAct->setFacetNormal(ivtmp);
321	delete ivtmp;
322	delete fAct;
323	}
324
325	/** \brief Default destructor
326	*/
327	~gcone()
328	{
329	//NOTE SAVE THE FACET STRUCTURE!!!
330	}
331
332
333	/** \brief Set the interior point of a cone */
334	void setIntPoint(intvec *iv)
335	{
336	this->ivIntPt=ivCopy(iv);
337	}
338
339	/** \brief Return the interior point */
340	intvec *getIntPoint()
341	{
342	return this->ivIntPt;
343	}
344
345	void showIntPoint()
346	{
347	ivIntPt->show();
348	}
349
350	void showFacets(short codim=1)
351	{
352	facet *f;
353	switch(codim)
354	{
355	case 1:
356	f = this->facetPtr;
357	break;
358	case 2:
359	f = this->facetPtr->codim2Ptr;
360	break;
361	}
362
363	intvec *iv = new intvec(this->numVars);
364
365	while (f->next!=NULL)
366	{
367	iv = f->getFacetNormal();
368	iv->show();
369	f=f->next;
370	}
371	//delete iv;
372	}
373
374	/** \brief Set gcone::numFacets */
375	void setNumFacets()
376	{
377	}
378
379	/** \brief Get gcone::numFacets */
380	int getNumFacets()
381	{
382	return this->numFacets;
383	}
384
385	int getUCN()
386	{
387	return this->UCN;
388	}
389
390	/** \brief Compute the normals of the cone
391	*
392	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
393	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
394	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
395	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
396	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
397	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
398	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
399	*
400	* Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
401	* an interior point of the cone.
402	*/
403	void getConeNormals(ideal const &I, bool compIntPoint=FALSE)
404	{
405	#ifdef gfan_DEBUG
406	std::cout << "* Computing Inequalities... *" << std::endl;
407	#endif
408	//All variables go here - except ineq matrix and *v, see below
409	int lengthGB=IDELEMS(I); // # of polys in the groebner basis
410	int pCompCount; // # of terms in a poly
411	poly aktpoly;
412	int numvar = pVariables; // # of variables in a polynomial (or ring?)
413	int leadexp[numvar]; // dirty hack of exp.vects
414	int aktexp[numvar];
415	int cols,rows; // will contain the dimensions of the ineq matrix - deprecated by
416	dd_rowrange ddrows;
417	dd_colrange ddcols;
418	dd_rowset ddredrows; // # of redundant rows in ddineq
419	dd_rowset ddlinset; // the opposite
420	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
421	dd_NumberType ddnumb=dd_Integer; //Number type
422	dd_ErrorType dderr=dd_NoError; //
423	// End of var declaration
424	#ifdef gfan_DEBUG
425	cout << "The Groebner basis has " << lengthGB << " elements" << endl;
426	cout << "The current ring has " << numvar << " variables" << endl;
427	#endif
428	cols = numvar;
429
430	//Compute the # inequalities i.e. rows of the matrix
431	rows=0; //Initialization
432	for (int ii=0;ii<IDELEMS(I);ii++)
433	{
434	aktpoly=(poly)I->m[ii];
435	rows=rows+pLength(aktpoly)-1;
436	}
437	#ifdef gfan_DEBUG
438	cout << "rows=" << rows << endl;
439	cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
440	#endif
441	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
442	dd_set_global_constants();
443	ddrows=rows;
444	ddcols=cols;
445	dd_MatrixPtr ddineq; //Matrix to store the inequalities
446	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
447
448	// We loop through each g\in GB and compute the resulting inequalities
449	for (int i=0; i<IDELEMS(I); i++)
450	{
451	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
452	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
453	#ifdef gfan_DEBUG
454	cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
455	#endif
456
457	int v=(int )omAlloc((numvar+1)*sizeof(int));
458	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
459
460	//Store leadexp for aktpoly
461	for (int kk=0;kk<numvar;kk++)
462	{
463	leadexp[kk]=v[kk+1];
464	//Since we need to know the difference of leadexp with the other expvects we do nothing here
465	//but compute the diff below
466	}
467
468
469	while (pNext(aktpoly)!=NULL) //move to next term until NULL
470	{
471	aktpoly=pNext(aktpoly);
472	pSetm(aktpoly); //doesn't seem to help anything
473	pGetExpV(aktpoly,v);
474	for (int kk=0;kk<numvar;kk++)
475	{
476	aktexp[kk]=v[kk+1];
477	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
478	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
479	}
480	aktmatrixrow=aktmatrixrow+1;
481	}//while
482
483	} //for
484
485	//Maybe add another row to contain the constraints of the standard simplex?
486
487	#ifdef gfan_DEBUG
488	cout << "The inequality matrix is" << endl;
489	dd_WriteMatrix(stdout, ddineq);
490	#endif
491
492	// The inequalities are now stored in ddineq
493	// Next we check for superflous rows
494	ddredrows = dd_RedundantRows(ddineq, &dderr);
495	if (dderr!=dd_NoError) // did an error occur?
496	{
497	dd_WriteErrorMessages(stderr,dderr); //if so tell us
498	} else
499	{
500	cout << "Redundant rows: ";
501	set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
502	}//if dd_Error
503
504	//Remove reduntant rows here!
505	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
506	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
507	ddcols = ddineq->colsize;
508	#ifdef gfan_DEBUG
509	cout << "Having removed redundancies, the normals now read:" << endl;
510	dd_WriteMatrix(stdout,ddineq);
511	cout << "Rows = " << ddrows << endl;
512	cout << "Cols = " << ddcols << endl;
513	#endif
514
515	/Write the normals into class facet/
516	#ifdef gfan_DEBUG
517	cout << "Creating list of normals" << endl;
518	#endif
519	/The pointer fRoot should be the return value of this function*/
520	facet *fRoot = new facet(); //instantiate new facet
521	this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
522	facet *fAct; //instantiate pointer to active facet
523	fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress!
524	//std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
525	for (int kk = 0; kk<ddrows; kk++)
526	{
527	intvec *load = new intvec(this->numVars); //intvec to store a single facet normal that will then be stored via setFacetNormal
528	for (int jj = 1; jj <ddcols; jj++)
529	{
530	double foo;
531	foo = mpq_get_d(ddineq->matrix[kk][jj]);
532	#ifdef gfan_DEBUG
533	std::cout << "fAct is " << foo << " at " << fAct << std::endl;
534	#endif
535	(*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
536	}//for (int jj = 1; jj <ddcols; jj++)
537
538	/Quick'n'dirty hack for flippability/
539	bool isFlippable=FALSE;
540	for (int jj = 0; jj<load->length(); jj++)
541	{
542	intvec *ivCanonical = new intvec(load->length());
543	(*ivCanonical)[jj]=1;
544	cout << "dotProd=" << dotProduct(load,ivCanonical) << endl;
545	if (dotProduct(load,ivCanonical)<0)
546	//if (ivMult(load,ivCanonical)<0)
547	{
548	isFlippable=TRUE;
549	break; //URGHS
550	}
551	}
552	if (isFlippable==FALSE)
553	{
554	cout << "Ignoring facet";
555	load->show();
556	//fAct->next=NULL;
557	}
558	else
559	{ /Now load should be full and we can call setFacetNormal/
560	fAct->setFacetNormal(load);
561	fAct->next = new facet();
562	fAct=fAct->next; //this should definitely not be called in the above while-loop :D
563	this->numFacets++;
564	}//if (isFlippable==FALSE)
565	//delete load;
566	}//for (int kk = 0; kk<ddrows; kk++)
567
568	/*
569	Now we should have a linked list containing the facet normals of those facets that are
570	-irredundant
571	-flipable
572	Adressing is done via *facetPtr
573	*/
574
575	if (compIntPoint==TRUE)
576	{
577	intvec *iv = new intvec(this->numVars);
578	interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets
579	this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt
580	//delete iv;
581	}
582
583	//Compute the number of facets
584	// wrong because ddineq->rowsize contains only irredundant facets. There can still be
585	// non-flippable ones!
586	//this->numFacets=ddineq->rowsize;
587
588	//Clean up but don't delete the return value! (Whatever it will turn out to be)
589	//dd_FreeMatrix(ddineq);
590	//set_free(ddredrows);
591	//free(ddnewpos);
592	//set_free(ddlinset);
593	//NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
594	//THIS SUCKS
595	//dd_free_global_constants();
596
597	}//method getConeNormals(ideal I)
598
599	/** \brief Compute the (codim-2)-facets of a given cone
600	* This method is used during noRevS
601	*/
602	void getCodim2Normals(gcone const &gc)
603	{
604	this->facetPtr->codim2Ptr = new facet(2); //instantiate a (codim-2)-facet
605	facet *codim2Act;
606	codim2Act = this->facetPtr->codim2Ptr;
607
608	dd_set_global_constants();
609
610	dd_MatrixPtr ddineq,P,ddakt;
611	dd_rowset impl_linset, redset;
612	dd_ErrorType err;
613	dd_rowindex newpos;
614
615	ddineq = facets2Matrix(gc); //get a matrix representation of the cone
616
617	/Now set appropriate linearity/
618	dd_PolyhedraPtr ddpolyh;
619	for (int ii=0; ii<this->numFacets; ii++)
620	{
621	ddakt = dd_CopyMatrix(ddineq);
622	set_addelem(ddakt->linset,ii+1);
623
624	dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
625
626	//dd_WriteMatrix(stdout,ddakt);
627	ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
628	P=dd_CopyGenerators(ddpolyh);
629	dd_WriteMatrix(stdout,P);
630
631	/* We loop through each row of P
632	* normalize it by making all entries integer ones
633	* and add the resulting vector to the int matrix facet::codim2Facets
634	*/
635	for (int jj=1;jj<=P->rowsize;jj++)
636	{
637	intvec *n = new intvec(this->numVars);
638	makeInt(P,jj,*n);
639	codim2Act->setFacetNormal(n);
640	codim2Act->next = new facet(2);
641	codim2Act = codim2Act->next;
642	delete n;
643	}
644
645	dd_FreeMatrix(ddakt);
646	dd_FreePolyhedra(ddpolyh);
647	}
648	}
649
650	/** \brief Compute the Groebner Basis on the other side of a shared facet
651	*
652	* Implements algorithm 4.3.2 from Anders' thesis.
653	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
654	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
655	* is parallel to \f$ leadexp(g) \f$
656	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
657	* Other possibilities include computing the rank of the matrix consisting of the vectors in question and
658	* computing an interior point of the facet and taking all terms having the same weight with respect
659	* to this interior point.
660	*\param ideal, facet
661	* Input: a marked,reduced Groebner basis and a facet
662	*/
663	void flip(ideal gb, facet *f) //Compute "the other side"
664	{
665	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
666	fNormal = f->getFacetNormal(); //read this->fNormal;
667	#ifdef gfan_DEBUG
668	std::cout << "===" << std::endl;
669	std::cout << "running gcone::flip" << std::endl;
670	// std::cout << "fNormal=";
671	// fNormal->show();
672	// std::cout << std::endl;
673	#endif
674	/1st step: Compute the initial ideal/
675	poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form
676	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
677	poly aktpoly;
678	intvec *check = new intvec(this->numVars); //array to store the difference of LE and v
679
680	for (int ii=0;ii<IDELEMS(gb);ii++)
681	{
682	aktpoly = (poly)gb->m[ii];
683	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
684	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
685	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
686	initialFormElement[ii]=pHead(aktpoly);
687
688	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
689	{
690	aktpoly=pNext(aktpoly); //next term
691	pSetm(aktpoly);
692	pGetExpV(aktpoly,v);
693	/* Convert (int)v into (intvec)check */
694	for (int jj=0;jj<this->numVars;jj++)
695	{
696	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
697	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
698	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
699	}
700	#ifdef gfan_DEBUG
701	// cout << "check=";
702	// check->show();
703	// cout << endl;
704	#endif
705	//TODO why not check, fNormal????
706	if (isParallel(check,fNormal)) //pass *check when
707	{
708	// cout << "Parallel vector found, adding to initialFormElement" << endl;
709	initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
710	}
711	}//while
712	#ifdef gfan_DEBUG
713	cout << "Initial Form=";
714	pWrite(initialFormElement[ii]);
715	cout << "---" << endl;
716	#endif
717	/Now initialFormElement must be added to (ideal)initialForm /
718	initialForm->m[ii]=initialFormElement[ii];
719	}//for
720	#ifdef gfan_DEBUG
721	cout << "Initial ideal is: " << endl;
722	idShow(initialForm);
723	//f->printFlipGB();
724	cout << "===" << endl;
725	#endif
726	//delete check;
727
728	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
729	/*Substep 2.1
730	compute $G_{-\alpha}(in_v(I))
731	see journal p. 66
732	*/
733	ring srcRing=currRing;
734
735	//intvec *negfNormal = new intvec(this->numVars);
736	//negfNormal=ivNeg(fNormal);
737	ring tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
738	rChangeCurrRing(tmpRing);
739
740	//rWrite(currRing); cout << endl;
741
742	ideal ina;
743	ina=idrCopyR(initialForm,srcRing);
744	#ifdef gfan_DEBUG
745	cout << "ina=";
746	idShow(ina); cout << endl;
747	#endif
748	ideal H;
749	H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
750	idSkipZeroes(H);
751	#ifdef gfan_DEBUG
752	// cout << "H="; idShow(H); cout << endl;
753	#endif
754	/*Substep 2.2
755	do the lifting and mark according to H
756	*/
757	rChangeCurrRing(srcRing);
758	ideal srcRing_H;
759	ideal srcRing_HH;
760	srcRing_H=idrCopyR(H,tmpRing);
761	#ifdef gfan_DEBUG
762	// cout << "srcRing_H = ";
763	// idShow(srcRing_H); cout << endl;
764	#endif
765	srcRing_HH=ffG(srcRing_H,this->gcBasis);
766	#ifdef gfan_DEBUG
767	// cout << "srcRing_HH = ";
768	// idShow(srcRing_HH); cout << endl;
769	#endif
770	/*Substep 2.2.1
771	Mark according to G_-\alpha
772	Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
773	we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
774	represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
775	Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
776	compute the difference accordingly
777	*/
778	dd_set_global_constants();
779	bool markingsAreCorrect=FALSE;
780	dd_MatrixPtr intPointMatrix;
781	int iPMatrixRows=0;
782	dd_rowrange aktrow=0;
783	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
784	{
785	poly aktpoly=(poly)srcRing_HH->m[ii];
786	iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
787	}
788	/* additionally one row for the standard-simplex and another for a row that becomes 0 during
789	construction of the differences
790	*/
791	intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
792	intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
793
794	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
795	{
796	markingsAreCorrect=FALSE; //crucial to initialise here
797	poly aktpoly=srcRing_HH->m[ii];
798	/Comparison of leading monomials is done via exponent vectors/
799	for (int jj=0;jj<IDELEMS(H);jj++)
800	{
801	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
802	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
803	pGetExpV(aktpoly,src_ExpV);
804	rChangeCurrRing(tmpRing); //this ring change is crucial!
805	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
806	rChangeCurrRing(srcRing);
807	bool expVAreEqual=TRUE;
808	for (int kk=1;kk<=this->numVars;kk++)
809	{
810	#ifdef gfan_DEBUG
811	//cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
812	#endif
813	if (src_ExpV[kk]!=dst_ExpV[kk])
814	{
815	expVAreEqual=FALSE;
816	}
817	}
818	//if (src_ExpV == dst_ExpV)
819	if (expVAreEqual==TRUE)
820	{
821	markingsAreCorrect=TRUE; //everything is fine
822	#ifdef gfan_DEBUG
823	// cout << "correct markings" << endl;
824	#endif
825	}//if (pHead(aktpoly)==pHead(H->m[jj])
826	delete src_ExpV;
827	delete dst_ExpV;
828	}//for (int jj=0;jj<IDELEMS(H);jj++)
829
830	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
831	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
832	if (markingsAreCorrect==TRUE)
833	{
834	pGetExpV(aktpoly,leadExpV);
835	}
836	else
837	{
838	rChangeCurrRing(tmpRing);
839	pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
840	rChangeCurrRing(srcRing);
841	}
842	/compute differences of the expvects/
843	while (pNext(aktpoly)!=NULL)
844	{
845	/*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
846	is not omitted when computing the differences*/
847	if(markingsAreCorrect==TRUE)
848	{
849	aktpoly=pNext(aktpoly);
850	pGetExpV(aktpoly,v);
851	}
852	else
853	{
854	pGetExpV(pHead(aktpoly),v);
855	markingsAreCorrect=TRUE;
856	}
857
858	for (int jj=0;jj<this->numVars;jj++)
859	{
860	/Store into ddMatrix/
861	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
862	}
863	aktrow +=1;
864	}
865	delete v;
866	delete leadExpV;
867	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
868	/Now we add the constraint for the standard simplex/
869	/NOTE:Might actually work without the standard simplex/
870	dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
871	for (int jj=1;jj<=this->numVars;jj++)
872	{
873	dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
874	}
875	dd_WriteMatrix(stdout,intPointMatrix);
876	intvec *iv_weight = new intvec(this->numVars);
877	interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point
878	dd_FreeMatrix(intPointMatrix);
879	dd_free_global_constants();
880
881	/*Step 3
882	turn the minimal basis into a reduced one
883	*/
884	int i,j;
885	ring dstRing=rCopy0(srcRing);
886	i=rBlocks(srcRing);
887
888	dstRing->order=(int )omAlloc((i+1)sizeof(int));
889	for(j=i;j>0;j--)
890	{
891	dstRing->order[j]=srcRing->order[j-1];
892	dstRing->block0[j]=srcRing->block0[j-1];
893	dstRing->block1[j]=srcRing->block1[j-1];
894	if (srcRing->wvhdl[j-1] != NULL)
895	{
896	dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
897	}
898	}
899	dstRing->order[0]=ringorder_a;
900	dstRing->order[1]=ringorder_dp;
901	dstRing->order[2]=ringorder_C;
902	dstRing->wvhdl[0] =( int )omAlloc((iv_weight->length())sizeof(int));
903
904	for (int ii=0;ii<this->numVars;ii++)
905	{
906	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
907	}
908	rComplete(dstRing);
909
910	// NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
911	// Thus:
912	//ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
913	//cout << "PLING" << endl;
914	/*ring dstRing=rCopy0(srcRing);
915	for (int ii=0;ii<this->numVars;ii++)
916	{
917	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
918	}*/
919	rChangeCurrRing(dstRing);
920	#ifdef gfan_DEBUG
921	rWrite(dstRing); cout << endl;
922	#endif
923	ideal dstRing_I;
924	dstRing_I=idrCopyR(srcRing_HH,srcRing);
925	//validOpts<1>=TRUE;
926	#ifdef gfan_DEBUG
927	//idShow(dstRing_I);
928	#endif
929	BITSET save=test;
930	test\|=Sy_bit(OPT_REDSB);
931	test\|=Sy_bit(6); //OPT_DEBUG
932	dstRing_I=kStd(idrCopyR(this->inputIdeal,this->rootRing),NULL,testHomog,NULL);
933	kInterRed(dstRing_I);
934	idSkipZeroes(dstRing_I);
935	test=save;
936	/End of step 3 - reduction/
937
938	f->setFlipGB(dstRing_I);//store the flipped GB
939	#ifdef gfan_DEBUG
940	cout << "Flipped GB is: " << endl;
941	f->printFlipGB();
942	#endif
943	}//void flip(ideal gb, facet *f)
944
945	/** \brief Compute the remainder of a polynomial by a given ideal
946	*
947	* Compute \f$ f^{\mathcal{G}} \f$
948	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
949	* However, since we are only interested in the remainder, there is no need to
950	* compute the factors \f$ a_i \f$
951	*/
952	//NOTE: Should be replaced by kNF or kNF2
953	poly restOfDiv(poly const &f, ideal const &I)
954	{
955	cout << "Entering restOfDiv" << endl;
956	poly p=f;
957	//pWrite(p);
958	//poly r=kCreateZeroPoly(,currRing,currRing); //The 0-polynomial, hopefully
959	poly r=NULL; //The zero polynomial
960	int ii;
961	bool divOccured;
962
963	while (p!=NULL)
964	{
965	ii=1;
966	divOccured=FALSE;
967
968	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
969	{
970	if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
971	{
972	poly step1,step2,step3;
973	//cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
974	step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
975	//cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
976	step2 = ppMult_qq(step1, I->m[ii-1]);
977	step3 = pSub(pCopy(p), step2);
978	//p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
979	//pSetm(p);
980	pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
981	p=step3;
982	//pWrite(p);
983	divOccured=TRUE;
984	}
985	else
986	{
987	ii += 1;
988	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
989	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
990	if (divOccured==FALSE)
991	{
992	//cout << "TICK 5" << endl;
993	r=pAdd(pCopy(r),pHead(p));
994	pSetm(r);
995	pSort(r);
996	//cout << "r="; pWrite(r); cout << endl;
997
998	if (pLength(p)!=1)
999	{
1000	p=pSub(pCopy(p),pHead(p)); //Here it may occur that p=0 instead of p=NULL
1001	}
1002	else
1003	{
1004	p=NULL; //Hack to correct this situation
1005	}
1006	//cout << "p="; pWrite(p);
1007	}//if (divOccured==FALSE)
1008	}//while (p!=0)
1009	return r;
1010	}//poly restOfDiv(poly const &f, ideal const &I)
1011
1012	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
1013	*/
1014	//NOTE: use kNF or kNF2 instead of restOfDivision
1015	ideal ffG(ideal const &H, ideal const &G)
1016	{
1017	cout << "Entering ffG" << endl;
1018	int size=IDELEMS(H);
1019	ideal res=idInit(size,1);
1020	poly temp1, temp2, temp3; //polys to temporarily store values for pSub
1021	for (int ii=0;ii<size;ii++)
1022	{
1023	res->m[ii]=restOfDiv(H->m[ii],G);
1024	//res->m[ii]=kNF(H->m[ii],G);
1025	temp1=H->m[ii];
1026	temp2=res->m[ii];
1027	temp3=pSub(temp1, temp2);
1028	res->m[ii]=temp3;
1029	//res->m[ii]=pSub(temp1,temp2); //buggy
1030	//pSort(res->m[ii]);
1031	//pSetm(res->m[ii]);
1032	//cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
1033	}
1034	return res;
1035	}
1036
1037	/** \brief Compute a Groebner Basis
1038	*
1039	* Computes the Groebner basis and stores the result in gcone::gcBasis
1040	*\param ideal
1041	*\return void
1042	*/
1043	void getGB(ideal const &inputIdeal)
1044	{
1045	ideal gb;
1046	gb=kStd(inputIdeal,NULL,testHomog,NULL);
1047	idSkipZeroes(gb);
1048	this->gcBasis=gb; //write the GB into gcBasis
1049	}//void getGB
1050
1051	/** \brief The Generic Groebner Walk due to FJLT
1052	* Needed for computing the search facet
1053	*/
1054	ideal GenGrbWlk(ideal, ideal)
1055	{
1056	}//GGW
1057
1058	/** \brief Compute the negative of a given intvec
1059	*/
1060	intvec *ivNeg(intvec const &iv)
1061	{
1062	intvec *res = new intvec(this->numVars);
1063	for(int ii=0;ii<this->numVars;ii++)
1064	{
1065	res[ii] = -(int)iv[ii];
1066	}
1067	return res;
1068	}
1069
1070
1071	/** \brief Compute the dot product of two intvecs
1072	*
1073	*/
1074	int dotProduct(intvec const &iva, intvec const &ivb)
1075	{
1076	//intvec iva=a;
1077	//intvec ivb=b;
1078	int res=0;
1079	for (int i=0;i<this->numVars;i++)
1080	{
1081	res = res+(iva[i]*ivb[i]);
1082	}
1083	return res;
1084	}//int dotProduct
1085
1086	/** \brief Check whether two intvecs are parallel
1087	*
1088	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1089	*/
1090	bool isParallel(intvec const &a, intvec const &b)
1091	{
1092	int lhs,rhs;
1093	lhs=dotProduct(a,b)*dotProduct(a,b);
1094	rhs=dotProduct(a,a)*dotProduct(b,b);
1095	//cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1096	if (lhs==rhs)
1097	{
1098	return TRUE;
1099	}
1100	else
1101	{
1102	return FALSE;
1103	}
1104	}//bool isParallel
1105
1106	/** \brief Compute an interior point of a given cone
1107	* Result will be written into intvec iv
1108	*/
1109	void interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
1110	{
1111	dd_LPPtr lp,lpInt;
1112	dd_ErrorType err=dd_NoError;
1113	dd_LPSolverType solver=dd_DualSimplex;
1114	dd_LPSolutionPtr lpSol=NULL;
1115	dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior
1116	dd_rowindex ddnewpos;
1117	dd_NumberType numb;
1118	//M->representation=dd_Inequality;
1119	//M->objective-dd_LPMin; //Not sure whether this is needed
1120
1121	//NOTE: Make this n-dimensional!
1122	//dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1123
1124	//dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1125	//if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1126	//cout << "Tick 2" << endl;
1127	//dd_WriteMatrix(stdout,M);
1128
1129	lp=dd_Matrix2LP(M, &err);
1130	if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}
1131	if (lp==NULL){cout << "LP is NULL" << endl;}
1132	#ifdef gfan_DEBUG
1133	// dd_WriteLP(stdout,lp);
1134	#endif
1135
1136	lpInt=dd_MakeLPforInteriorFinding(lp);
1137	if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1138	#ifdef gfan_DEBUG
1139	// dd_WriteLP(stdout,lpInt);
1140	#endif
1141
1142	dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1143	if (err!=dd_NoError)
1144	{
1145	cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1146	dd_WriteErrorMessages(stdout, err);
1147	}
1148
1149	//dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior
1150	if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1151	//cout << "Tick 5" << endl;
1152
1153	//lpSol=dd_CopyLPSolution(lpInt);
1154	if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}
1155	//cout << "Tick 6" << endl;
1156	#ifdef gfan_DEBUG
1157	cout << "Interior point: ";
1158	#endif
1159	for (int ii=1; ii<(lpSol->d)-1;ii++)
1160	{
1161	#ifdef gfan_DEBUG
1162	dd_WriteNumber(stdout,lpSol->sol[ii]);
1163	#endif
1164	/* NOTE This works only as long as gmp returns fractions with the same denominator*/
1165	(iv)[ii-1]=(int)mpz_get_d(mpq_numref(lpSol->sol[ii])); //looks evil, but does the trick
1166	}
1167	dd_FreeLPSolution(lpSol);
1168	dd_FreeLPData(lpInt);
1169	dd_FreeLPData(lp);
1170	set_free(ddlinset);
1171	set_free(ddredrows);
1172
1173	}//void interiorPoint(dd_MatrixPtr const &M)
1174
1175	/** \brief Copy a ring and add a weighted ordering in first place
1176	* Kudos to walkSupport.cc
1177	*/
1178	ring rCopyAndAddWeight(ring const &r, intvec const *ivw)
1179	{
1180	ring res=(ring)omAllocBin(ip_sring_bin);
1181	memcpy4(res,r,sizeof(ip_sring));
1182	res->VarOffset = NULL;
1183	res->ref=0;
1184
1185	if (r->algring!=NULL)
1186	r->algring->ref++;
1187	if (r->parameter!=NULL)
1188	{
1189	res->minpoly=nCopy(r->minpoly);
1190	int l=rPar(r);
1191	res->parameter=(char *)omAlloc(lsizeof(char_ptr));
1192	int i;
1193	for(i=0;i<rPar(r);i++)
1194	{
1195	res->parameter[i]=omStrDup(r->parameter[i]);
1196	}
1197	}
1198
1199	int i=rBlocks(r);
1200	int jj;
1201
1202	res->order =(int )omAlloc((i+1)sizeof(int));
1203	res->block0=(int )omAlloc((i+1)sizeof(int));
1204	res->block1=(int )omAlloc((i+1)sizeof(int));
1205	res->wvhdl =(int *)omAlloc((i+1)sizeof(int**));
1206	for(jj=0;jj<i;jj++)
1207	{
1208	if (r->wvhdl[jj] != NULL)
1209	{
1210	res->wvhdl[jj] = (int*) omMemDup(r->wvhdl[jj-1]);
1211	}
1212	else
1213	{
1214	res->wvhdl[jj+1]=NULL;
1215	}
1216	}
1217
1218	for (jj=0;jj<i;jj++)
1219	{
1220	res->order[jj+1]=r->order[jj];
1221	res->block0[jj+1]=r->block0[jj];
1222	res->block1[jj+1]=r->block1[jj];
1223	}
1224
1225	res->order[0]=ringorder_a;
1226	res->order[1]=ringorder_dp; //basically useless, since that should never be used
1227	int length=ivw->length();
1228	int A=(int )omAlloc(length*sizeof(int));
1229	for (jj=0;jj<length;jj++)
1230	{
1231	A[jj]=(*ivw)[jj];
1232	}
1233	res->wvhdl[0]=(int *)A;
1234	res->block0[0]=1;
1235	res->block1[0]=length;
1236
1237	res->names = (char *)omAlloc0(rVar(r) sizeof(char_ptr));
1238	for (i=rVar(res)-1;i>=0; i--)
1239	{
1240	res->names[i] = omStrDup(r->names[i]);
1241	}
1242	rComplete(res);
1243	return res;
1244	}//rCopyAndAdd
1245
1246	ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1247	{
1248	ring res=rCopy0(currRing);
1249	rComplete(res);
1250	rSetWeightVec(res,(int64*)ivw);
1251	//rChangeCurrRing(rnew);
1252	return res;
1253	}
1254
1255	/** \brief Checks whether a given facet is a search facet
1256	* Determines whether a given facet of a cone is the search facet of a neighbouring cone
1257	* This is done in the following way:
1258	* We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1259	* that is first crossed during the generic walk.
1260	* We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1261	* If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1262	*/
1263	bool isSearchFacet(gcone &gcTmp, facet *testfacet)
1264	{
1265	ring actRing=currRing;
1266	facet facetPtr=(facet)gcTmp.facetPtr;
1267	facet fMin=new facet(facetPtr); //Pointer to the "minimal" facet
1268	//facet *fMin = new facet(tmpcone.facetPtr);
1269	//fMin=tmpcone.facetPtr; //Initialise to first facet of tmpcone
1270	facet *fAct; //Ptr to alpha_i
1271	facet *fCmp; //Ptr to alpha_j
1272	fAct = fMin;
1273	fCmp = fMin->next;
1274
1275	rChangeCurrRing(this->rootRing); //because we compare the monomials in the rootring
1276	poly p=pInit();
1277	poly q=pInit();
1278	intvec *alpha_i = new intvec(this->numVars);
1279	intvec *alpha_j = new intvec(this->numVars);
1280	intvec *sigma = new intvec(this->numVars);
1281	sigma=gcTmp.getIntPoint();
1282
1283	int u=(int )omAlloc((this->numVars+1)*sizeof(int));
1284	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1285	u[0]=0; v[0]=0;
1286	int weight1,weight2;
1287	while(fAct->next->next!=NULL) //NOTE this is ugly. Can it be done without fCmp?
1288	{
1289	/* Get alpha_i and alpha_{i+1} */
1290	alpha_i=fAct->getFacetNormal();
1291	alpha_j=fCmp->getFacetNormal();
1292	#ifdef gfan_DEBUG
1293	alpha_i->show();
1294	alpha_j->show();
1295	#endif
1296	/Compute the dot product of sigma and alpha_{i,j}/
1297	weight1=dotProduct(sigma,alpha_i);
1298	weight2=dotProduct(sigma,alpha_j);
1299	#ifdef gfan_DEBUG
1300	cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1301	#endif
1302	/Adjust alpha_i and alpha_i+1 accordingly/
1303	for(int ii=1;ii<=this->numVars;ii++)
1304	{
1305	u[ii]=weight1(alpha_i)[ii-1];
1306	v[ii]=weight2(alpha_j)[ii-1];
1307	}
1308
1309	/Now p_weight and q_weight need to be compared as exponent vectors/
1310	pSetCoeff0(p,nInit(1));
1311	pSetCoeff0(q,nInit(1));
1312	pSetExpV(p,u);
1313	pSetm(p);
1314	pSetExpV(q,v);
1315	pSetm(q);
1316	#ifdef gfan_DEBUG
1317	pWrite(p);pWrite(q);
1318	#endif
1319	/*We want to check whether x^p < x^q
1320	=> want to check for return value 1 */
1321	if (pLmCmp(p,q)==1) //i.e. x^q is smaller
1322	{
1323	fMin=fCmp;
1324	fAct=fMin;
1325	fCmp=fCmp->next;
1326	}
1327	else
1328	{
1329	//fAct=fAct->next;
1330	if(fCmp->next!=NULL)
1331	{
1332	fCmp=fCmp->next;
1333	}
1334	else
1335	{
1336	fAct=fAct->next;
1337	}
1338	}
1339	//fAct=fAct->next;
1340	}//while(fAct.facetPtr->next!=NULL)
1341	delete alpha_i,alpha_j,sigma;
1342
1343	/If testfacet was minimal then fMin should still point there /
1344
1345	//if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))
1346	#ifdef gfan_DEBUG
1347	cout << "Checking for parallelity" << endl <<" fMin is";
1348	fMin->printNormal();
1349	cout << "testfacet is ";
1350	testfacet->printNormal();
1351	cout << endl;
1352	#endif
1353	if (fMin==gcTmp.facetPtr)
1354	//if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
1355	//if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
1356	{
1357	cout << "Parallel" << endl;
1358	rChangeCurrRing(actRing);
1359	//delete alpha_min, test;
1360	return TRUE;
1361	}
1362	else
1363	{
1364	cout << "Not parallel" << endl;
1365	rChangeCurrRing(actRing);
1366	//delete alpha_min, test;
1367	return FALSE;
1368	}
1369	}//bool isSearchFacet
1370
1371	/** \brief Check for equality of two intvecs
1372	*/
1373	bool areEqual(intvec const &a, intvec const &b)
1374	{
1375	bool res=TRUE;
1376	for(int ii=0;ii<this->numVars;ii++)
1377	{
1378	if(a[ii]!=b[ii])
1379	{
1380	res=FALSE;
1381	break;
1382	}
1383	}
1384	return res;
1385	}
1386
1387	/** \brief The reverse search algorithm
1388	*/
1389	void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1390	{
1391	facet *fAct=new facet();
1392	fAct = gcAct->facetPtr;
1393
1394	while(fAct->next!=NULL) //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
1395	{
1396	cout << "==========================================================================================="<< endl;
1397	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1398	gcone gcTmp = new gcone(gcAct);
1399	//idShow(gcTmp->gcBasis);
1400	gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1401	#ifdef gfan_DEBUG
1402	facet *f = new facet();
1403	f=gcTmp->facetPtr;
1404	while(f->next!=NULL)
1405	{
1406	f->printNormal();
1407	f=f->next;
1408	}
1409	#endif
1410	gcTmp->showIntPoint();
1411	/recursive part goes gere/
1412	if (isSearchFacet(gcTmp,(facet)gcAct->facetPtr))
1413	{
1414	gcAct->next=gcTmp;
1415	cout << "PING"<< endl;
1416	reverseSearch(gcTmp);
1417	}
1418	else
1419	{
1420	delete gcTmp;
1421	/NOTE remove fAct from linked list. It's no longer needed/
1422	}
1423	/recursion ends/
1424	fAct = fAct->next;
1425	}//while(fAct->next!=NULL)
1426	}//reverseSearch
1427
1428	/** \brief The new method of Markwig and Jensen
1429	* Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1430	* In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1431	* Keep a list of facets as a linked list containing an intvec and an integer matrix.
1432	* Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1433	* soon as we compute this facet again. Comparison of facets is done by...
1434	*/
1435	void noRevS(gcone &gc, bool usingIntPoint=FALSE)
1436	{
1437	facet *fListPtr = new facet(); //The list containing ALL facets we come across
1438	facet *fAct;// = new facet();
1439	fAct = gc.facetPtr;
1440	fListPtr = gc.facetPtr;
1441	#ifdef gfan_DEBUG
1442	cout << "NoRevs" << endl;
1443	cout << "Facets are:" << endl;
1444	gc.showFacets();
1445	#endif
1446	//dd_set_global_constants();
1447	dd_rowrange ddrows;
1448	dd_colrange ddcols;
1449	ddrows=2; //Each facet is described by its normal
1450	ddcols=gc.numVars+1; // plus one for the standard simplex
1451	if(usingIntPoint){
1452	while(fAct->next!=NULL)
1453	{
1454	/Compute an interior point for each facet/
1455	dd_MatrixPtr ddineq;
1456	ddineq=dd_CreateMatrix(ddrows,ddcols);
1457	intvec *comp = new intvec(this->numVars);
1458	comp=fAct->getFacetNormal();
1459	for(int ii=0; ii<this->numVars; ii++)
1460	{
1461	dd_set_si(ddineq->matrix[0][ii+1],(*comp)[ii]);
1462	dd_set_si(ddineq->matrix[1][ii+1],1); //Assure we search in the pos. orthant
1463	}
1464	set_addelem(ddineq->linset,1); //We want equality in the first row
1465	//dd_WriteMatrix(stdout,ddineq);
1466	interiorPoint(ddineq,*comp);
1467	/**/
1468	#ifdef gfan_DEBUG
1469	cout << "IP is";
1470	comp->show(); cout << endl;
1471	#endif
1472	//Store the interior point and the UCN
1473	fListPtr->setInteriorPoint( comp );
1474	fListPtr->setUCN( gc.getUCN() );
1475
1476	dd_FreeMatrix(ddineq);
1477	fAct=fAct->next; //iterate
1478	}
1479	}
1480	else/Facets of facets/
1481	{
1482
1483	this->getCodim2Normals(gc);
1484
1485	//Compute unique representation of codim-2-facets
1486	this->normalize();
1487
1488	gc.writeConeToFile(gc);
1489
1490	/*2nd step
1491	Choose a facet from fListPtr, flip it and forget the previous cone
1492	We always choose the first facet from fListPtr as facet to be flipped
1493	*/
1494
1495
1496
1497	}//else usingIntPoint
1498
1499
1500
1501
1502	//NOTE Hm, comment in and get a crash for free...
1503	//dd_free_global_constants();
1504	//gc.writeConeToFile(gc);
1505
1506
1507	//fAct = fListPtr;
1508	//gcone *gcTmp = new gcone(gc); //copy
1509	//gcTmp->flip(gcTmp->gcBasis,fAct);
1510	//NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
1511
1512	}//void noRevS(gcone &gc)
1513
1514
1515	/** \brief Make a set of rational vectors into integers
1516	*
1517	* We compute the lcm of the denominators and multiply with this to get integer values.
1518	* \param dd_MatrixPtr,intvec
1519	*/
1520	void makeInt(dd_MatrixPtr const &M, int const line, intvec &n)
1521	{
1522	mpz_t denom[this->numVars];
1523	for(int ii=0;ii<this->numVars;ii++)
1524	{
1525	mpz_init_set_str(denom[ii],"0",10);
1526	}
1527
1528	mpz_t kgV,tmp;
1529	mpz_init(kgV);
1530	mpz_init(tmp);
1531
1532	for (int ii=0;ii<(M->colsize)-1;ii++)
1533	{
1534	mpz_t z;
1535	mpz_init(z);
1536	mpq_get_den(z,M->matrix[line-1][ii+1]);
1537	mpz_set( denom[ii], z);
1538	mpz_clear(z);
1539	}
1540
1541	/Compute lcm of the denominators/
1542	mpz_set(tmp,denom[0]);
1543	for (int ii=0;ii<(M->colsize)-1;ii++)
1544	{
1545	mpz_lcm(kgV,tmp,denom[ii]);
1546	}
1547
1548	/Multiply the nominators by kgV/
1549	mpq_t qkgV,res;
1550	mpq_init(qkgV);
1551	mpq_set_str(qkgV,"1",10);
1552	mpq_canonicalize(qkgV);
1553
1554	mpq_init(res);
1555	mpq_set_str(res,"1",10);
1556	mpq_canonicalize(res);
1557
1558	mpq_set_num(qkgV,kgV);
1559
1560	// mpq_canonicalize(qkgV);
1561	for (int ii=0;ii<(M->colsize)-1;ii++)
1562	{
1563	mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
1564	n[ii]=(int)mpz_get_d(mpq_numref(res));
1565	}
1566	//mpz_clear(denom[this->numVars]);
1567	mpz_clear(kgV);
1568	mpq_clear(qkgV); mpq_clear(res);
1569
1570	}
1571	/**
1572	* We compute the gcd of the components of the codim-2-facets and
1573	* multiply the each codim-2facet by it. Thus we get a normalized representation of each
1574	* (codim-2)-facet normal.
1575	*/
1576	void normalize()
1577	{
1578	int ggT=1;
1579	facet *codim2Act;
1580	codim2Act = this->facetPtr->codim2Ptr;
1581	intvec *n = new intvec(this->numVars);
1582
1583	while(codim2Act->next!=NULL)
1584	{
1585	n=codim2Act->getFacetNormal();
1586	for(int ii=0;ii<this->numVars;ii++)
1587	{
1588	ggT = intgcd(ggT,(*n)[ii]);
1589	}
1590	codim2Act = codim2Act->next;
1591	}
1592	//delete n;
1593	codim2Act = this->facetPtr->codim2Ptr; //reset to start of linked list
1594	while(codim2Act->next!=NULL)
1595	{
1596	//intvec *n = new intvec(this->numVars);
1597	n=codim2Act->getFacetNormal();
1598	intvec *new_n = new intvec(this->numVars);
1599	for(int ii=0;ii<this->numVars;ii++)
1600	{
1601	(new_n)[ii] = (n)[ii]*ggT;
1602	}
1603	codim2Act->setFacetNormal(new_n);
1604	codim2Act = codim2Act->next;
1605	//delete n;
1606	//delete new_n;
1607	}
1608	}
1609
1610	/** \brief Compute the gcd of two ints
1611	*/
1612	int intgcd(int a, int b)
1613	{
1614	int r, p=a, q=b;
1615	if(p < 0)
1616	p = -p;
1617	if(q < 0)
1618	q = -q;
1619	while(q != 0)
1620	{
1621	r = p % q;
1622	p = q;
1623	q = r;
1624	}
1625	return p;
1626	}
1627
1628	/** \brief Construct a dd_MatrixPtr from a cone's list of facets
1629	*
1630	*/
1631	dd_MatrixPtr facets2Matrix(gcone const &gc)
1632	{
1633	facet *fAct;
1634	fAct = gc.facetPtr;
1635
1636	dd_MatrixPtr M;
1637	dd_rowrange ddrows;
1638	dd_colrange ddcols;
1639	ddcols=(this->numVars)+1;
1640	ddrows=this->numFacets;
1641	dd_NumberType numb = dd_Integer;
1642	M=dd_CreateMatrix(ddrows,ddcols);
1643
1644	//dd_set_global_constants();
1645
1646	//intvec *comp = new intvec(this->numVars);
1647
1648	/*for (int jj=0; jj<M->rowsize; jj++)
1649	{
1650	intvec *comp = new intvec(this->numVars);
1651	comp = fAct->getFacetNormal();
1652	for (int ii=0; ii<this->numVars; ii++)
1653	{
1654	dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
1655	}
1656	if(fAct->next!=NULL)
1657	{
1658	fAct = fAct->next;
1659	}
1660	delete comp;
1661	}//jj
1662	*/
1663	int jj=0;
1664	while(fAct->next!=NULL)
1665	{
1666	intvec *comp;
1667	comp = fAct->getFacetNormal();
1668	for(int ii=0;ii<this->numVars;ii++)
1669	{
1670	dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
1671	}
1672	jj++;
1673	fAct=fAct->next;
1674	}
1675
1676	return M;
1677	}
1678
1679	/** \brief Write information about a cone into a file on disk
1680	*
1681	* This methods writes the information needed for the "second" method into a file.
1682	* The file's is divided in sections containing information on
1683	* 1) the ring
1684	* 2) the cone's Groebner Basis
1685	* 3) the cone's facets
1686	* Each line contains exactly one date
1687	* Each section starts with its name in CAPITALS
1688	*/
1689	void writeConeToFile(gcone const &gc, bool usingIntPoints=FALSE)
1690	{
1691	ofstream gcOutputFile("/tmp/cone1.gc");
1692	facet *fAct = new facet();
1693	fAct = gc.facetPtr;
1694
1695	if (!gcOutputFile)
1696	{
1697	cout << "Error opening file for writing in writeConeToFile" << endl;
1698	}
1699	else
1700	{ /*gcOutputFile << "UCN" << endl;
1701	gcOutputFile << this->UCN << endl;*/
1702	gcOutputFile << "RING" << endl;
1703	//Write this->gcBasis into file
1704	gcOutputFile << "GCBASIS" << endl;
1705	for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
1706	{
1707	gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
1708	}
1709
1710	gcOutputFile << "FACETS" << endl;
1711	while(fAct->next!=NULL)
1712	{
1713	intvec *iv = new intvec(gc.numVars);
1714	iv=fAct->getFacetNormal();
1715	for (int ii=0;ii<iv->length();ii++)
1716	{
1717	if (ii<iv->length()-1)
1718	{
1719	gcOutputFile << (*iv)[ii] << ",";
1720	}
1721	else
1722	{
1723	gcOutputFile << (*iv)[ii] << endl;
1724	}
1725	}
1726	fAct=fAct->next;
1727	//delete iv; iv=NULL;
1728	}
1729
1730	gcOutputFile.close();
1731	//delete fAct; fAct=NULL;
1732	}
1733
1734	}//writeConeToFile(gcone const &gc)
1735
1736	/** \brief Reads a cone from a file identified by its number
1737	*/
1738	void readConeFromFile(int gcNum)
1739	{
1740	}
1741
1742	friend class facet;
1743	};//class gcone
1744
1745	int gcone::UCN=0;
1746
1747	ideal gfan(ideal inputIdeal)
1748	{
1749	int numvar = pVariables;
1750
1751	enum searchMethod {
1752	reverseSearch,
1753	noRevS
1754	};
1755
1756	searchMethod method;
1757	method = noRevS;
1758	// method = reverseSearch;
1759
1760	#ifdef gfan_DEBUG
1761	cout << "Now in subroutine gfan" << endl;
1762	#endif
1763	ring inputRing=currRing; // The ring the user entered
1764	ring rootRing; // The ring associated to the target ordering
1765	ideal res;
1766	facet *fRoot;
1767
1768	if (method==reverseSearch)
1769	{
1770
1771	/* Construct a new ring which will serve as our root*/
1772	rootRing=rCopy0(currRing);
1773	rootRing->order[0]=ringorder_lp;
1774
1775	rComplete(rootRing);
1776	rChangeCurrRing(rootRing);
1777
1778	/* Fetch the inputIdeal into our rootRing */
1779	map theMap=(map)idMaxIdeal(1); //evil hack!
1780	theMap->preimage=NULL; //neccessary?
1781	ideal rootIdeal;
1782	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
1783	#ifdef gfan_DEBUG
1784	cout << "Root ideal is " << endl;
1785	idShow(rootIdeal);
1786	cout << "The root ring is " << endl;
1787	rWrite(rootRing);
1788	cout << endl;
1789	#endif
1790
1791	//gcone *gcRoot = new gcone(); //Instantiate the sink
1792	gcone *gcRoot = new gcone(rootRing,rootIdeal);
1793	gcone *gcAct;
1794	gcAct = gcRoot;
1795	gcAct->numVars=pVariables;
1796	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
1797	idShow(gcAct->gcBasis);
1798	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
1799	//gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1800	/*Now it is time to compute the search facets, respectively start the reverse search.
1801	But since we are in the root all facets should be search facets. IS THIS TRUE?
1802	NOTE: Check for flippability is not very sophisticated
1803	*/
1804	//gcAct->reverseSearch(gcAct);
1805	rChangeCurrRing(rootRing);
1806	res=gcRoot->gcBasis;
1807	}//if method==reverSearch
1808
1809	if(method==noRevS)
1810	{
1811	gcone *gcRoot = new gcone(currRing,inputIdeal);
1812	gcone *gcAct;
1813	gcAct = gcRoot;
1814	gcAct->numVars=pVariables;
1815	gcAct->getGB(inputIdeal);
1816	gcAct->getConeNormals(gcAct->gcBasis);
1817	gcAct->noRevS(*gcAct);
1818	res=gcAct->gcBasis;
1819	}
1820
1821	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
1822	The return type will then be of type LIST_CMD
1823	Assume gfan has finished, thus we have enumerated all the cones
1824	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
1825	Groebner Basis and merge this somehow with LIST_CMD
1826	=> Count the cones!
1827	*/
1828	//rChangeCurrRing(rootRing);
1829	//res=gcAct->gcBasis;
1830	//res=gcRoot->gcBasis;
1831	return res;
1832	//return GBlist;
1833	}
1834	/*
1835	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
1836	*/
1837	#endif

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

Download in other formats: