Context Navigation

source: git/kernel/gfan.cc @ 52696b

spielwiese

Last change on this file since 52696b was 52696b, checked in by Martin Monerjan, 14 years ago
minor fix in normalize git-svn-id: file:///usr/local/Singular/svn/trunk@11942 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 61.2 KB

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

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

Download in other formats: