Context Navigation

source: git/kernel/gfan.cc @ d5042e

spielwiese

Last change on this file since d5042e was d5042e, checked in by Martin Monerjan, 14 years ago
Modified facet::facet() to set everything to 0 in the beginning New method enqueueNewFacets BUG in the way getConeNormals creates the list of normals. There is one facet to much created everytime! git-svn-id: file:///usr/local/Singular/svn/trunk@11915 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 58.3 KB

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

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

Download in other formats: