Context Navigation

source: git/kernel/gfan.cc @ 6e51338

spielwiese

Last change on this file since 6e51338 was 6e51338, checked in by Martin Monerjan, 14 years ago
Keep matrix of facet normals in each gcone in order to have the right codim2facets git-svn-id: file:///usr/local/Singular/svn/trunk@11943 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 61.7 KB

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

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

Download in other formats: