Context Navigation

source: git/kernel/gfan.cc @ b7510d

spielwiese

Last change on this file since b7510d was b7510d, checked in by Martin Monerjan, 14 years ago
BUGFIX: Creation of linked list of facets now working correctly New method showSLA() - debugging purposes only Fixed enumeration of cones via UCN -> works now Enqueueing into SLA done TODO: Find out why it crashes in flip(). Maybe the usual wrong ring issue git-svn-id: file:///usr/local/Singular/svn/trunk@11920 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 60.2 KB

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

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

Download in other formats: