Context Navigation

source: git/kernel/gfan.cc @ 005d00a

spielwiese

Last change on this file since 005d00a was 005d00a, checked in by Martin Monerjan, 14 years ago
have return type LIST_CMD git-svn-id: file:///usr/local/Singular/svn/trunk@12283 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 81.6 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009/11/03 06:57:32 $
5	$Header: /usr/local/Singular/cvsroot/kernel/gfan.cc,v 1.103 2009/11/03 06:57:32 monerjan Exp $
6	$Id$
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 "structs.h"
23	#include "../Singular/lists.h"
24	#include "prCopy.h"
25	#include <iostream>
26	#include <bitset>
27	#include <fstream> //read-write cones to files
28	#include <gmp.h>
29	#include <string>
30	#include <sstream>
31	#include <time.h>
32
33	//#include <gmpxx.h>
34
35	/DO NOT REMOVE THIS/
36	#ifndef GMPRATIONAL
37	#define GMPRATIONAL
38	#endif
39
40	//Hacks for different working places
41	#define p800
42
43	#ifdef UNI
44	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
45	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
46	#endif
47
48	#ifdef HOME
49	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
50	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
51	#endif
52
53	#ifdef ITWM
54	#include "/u/slg/monerjan/cddlib/include/setoper.h"
55	#include "/u/slg/monerjan/cddlib/include/cdd.h"
56	#include "/u/slg/monerjan/cddlib/include/cddmp.h"
57	#endif
58
59	#ifdef p800
60	#include "../../cddlib/include/setoper.h"
61	#include "../../cddlib/include/cdd.h"
62	#include "../../cddlib/include/cddmp.h"
63	#endif
64
65	#ifndef gfan_DEBUG
66	//#define gfan_DEBUG
67	#ifndef gfan_DEBUGLEVEL
68	#define gfan_DEBUGLEVEL 1
69	#endif
70	#endif
71
72	#include <gfan.h>
73	using namespace std;
74
75	/**
76	*\brief Class facet
77	* Implements the facet structure as a linked list
78	*
79	*/
80
81	/** \brief The default constructor for facets
82	* Do I need a constructor of type facet(intvec)?
83	*/
84	facet::facet()
85	{
86	// Pointer to next facet. */
87	/* Defaults to NULL. This way there is no need to check explicitly */
88	this->fNormal=NULL;
89	this->interiorPoint=NULL;
90	this->UCN=0;
91	this->codim2Ptr=NULL;
92	this->codim=1; //default for (codim-1)-facets
93	this->numCodim2Facets=0;
94	this->flipGB=NULL;
95	this->isIncoming=FALSE;
96	this->next=NULL;
97	this->prev=NULL;
98	this->flipRing=NULL; //the ring on the other side
99	this->isFlippable=FALSE;
100	}
101
102	/** \brief Constructor for facets of codim >= 2
103	* Note that as of now the code of the constructors is only for facets and codim2-faces. One
104	* could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
105	*/
106	facet::facet(int const &n)
107	{
108	this->fNormal=NULL;
109	this->interiorPoint=NULL;
110	this->UCN=0;
111	this->codim2Ptr=NULL;
112	if(n>1)
113	{
114	this->codim=n;
115	}//NOTE Handle exception here!
116	this->numCodim2Facets=0;
117	this->flipGB=NULL;
118	this->isIncoming=FALSE;
119	this->next=NULL;
120	this->prev=NULL;
121	this->flipRing=NULL;
122	this->isFlippable=FALSE;
123	}
124
125	/** \brief The copy constructor
126	*/
127	facet::facet(const facet& f)
128	{
129	}
130
131	/** The default destructor */
132	facet::~facet()
133	{
134	delete this->fNormal;
135	delete this->interiorPoint;
136	/* Cleanup the codim2-structure */
137	if(this->codim==2)
138	{
139	facet *codim2Ptr;
140	codim2Ptr = this->codim2Ptr;
141	while(codim2Ptr!=NULL)
142	{
143	delete codim2Ptr->fNormal;
144	codim2Ptr = codim2Ptr->next;
145	}
146	}
147	delete this->codim2Ptr;
148	idDelete((ideal *)&this->flipGB);
149	rDelete(this->flipRing);
150	this->flipRing=NULL;
151	//this=NULL;
152	}
153
154
155	/** \brief Comparison of facets*/
156	bool facet::areEqual(facet &f, facet &g)
157	{
158	bool res = TRUE;
159	intvec *fNormal = new intvec(pVariables);
160	intvec *gNormal = new intvec(pVariables);
161	fNormal = f.getFacetNormal();
162	gNormal = g.getFacetNormal();
163	if((fNormal == gNormal))//\|\|(gcone::isParallel(fNormal,gNormal)))
164	{
165	if(f.numCodim2Facets==g.numCodim2Facets)
166	{
167	facet *f2Act;
168	facet *g2Act;
169	f2Act = f.codim2Ptr;
170	g2Act = g.codim2Ptr;
171	intvec *f2N = new intvec(pVariables);
172	intvec *g2N = new intvec(pVariables);
173	while(f2Act->next!=NULL && g2Act->next!=NULL)
174	{
175	for(int ii=0;ii<f2N->length();ii++)
176	{
177	if(f2Act->getFacetNormal() != g2Act->getFacetNormal())
178	{
179	res=FALSE;
180	}
181	if (res==FALSE)
182	break;
183	}
184	if(res==FALSE)
185	break;
186
187	f2Act = f2Act->next;
188	g2Act = g2Act->next;
189	}//while
190	}//if
191	}//if
192	else
193	{
194	res = FALSE;
195	}
196	return res;
197	}
198
199	/** Stores the facet normal \param intvec*/
200	void facet::setFacetNormal(intvec *iv)
201	{
202	this->fNormal = ivCopy(iv);
203	}
204
205	/** Hopefully returns the facet normal */
206	intvec *facet::getFacetNormal()
207	{
208	return this->fNormal;
209	}
210
211	/** Method to print the facet normal*/
212	void facet::printNormal()
213	{
214	fNormal->show();
215	}
216
217	/** Store the flipped GB*/
218	void facet::setFlipGB(ideal I)
219	{
220	this->flipGB=idCopy(I);
221	}
222
223	/** Return the flipped GB*/
224	ideal facet::getFlipGB()
225	{
226	return this->flipGB;
227	}
228
229	/** Print the flipped GB*/
230	void facet::printFlipGB()
231	{
232	#ifndef NDEBUG
233	idShow(this->flipGB);
234	#endif
235	}
236
237	/** Set the UCN */
238	void facet::setUCN(int n)
239	{
240	this->UCN=n;
241	}
242
243	/** \brief Get the UCN
244	* Returns the UCN iff this != NULL, else -1
245	*/
246	int facet::getUCN()
247	{
248	if(this!=NULL && ( this!=(facet )0xfbfbfbfb \|\| this!=(facet )0xfbfbfbfbfbfbfbfb) )
249	return this->UCN;
250	else
251	return -1;
252	}
253
254	/** Store an interior point of the facet */
255	void facet::setInteriorPoint(intvec *iv)
256	{
257	this->interiorPoint = ivCopy(iv);
258	}
259
260	intvec *facet::getInteriorPoint()
261	{
262	return this->interiorPoint;
263	}
264
265	/** \brief Debugging function
266	* prints the facet normal an all (codim-2)-facets that belong to it
267	*/
268	void facet::fDebugPrint()
269	{
270	facet *codim2Act;
271	codim2Act = this->codim2Ptr;
272	intvec *fNormal;
273	fNormal = this->getFacetNormal();
274	intvec *f2Normal;
275	cout << "=======================" << endl;
276	cout << "Facet normal = (";
277	fNormal->show(1,1);
278	cout << ")"<<endl;
279	cout << "-----------------------" << endl;
280	cout << "Codim2 facets:" << endl;
281	while(codim2Act!=NULL)
282	{
283	f2Normal = codim2Act->getFacetNormal();
284	cout << "(";
285	f2Normal->show(1,0);
286	cout << ")" << endl;
287	codim2Act = codim2Act->next;
288	}
289	cout << "=======================" << endl;
290	}
291
292	//friend class gcone; //Bad style
293
294
295	/**
296	*\brief Implements the cone structure
297	*
298	* A cone is represented by a linked list of facet normals
299	* @see facet
300	*/
301
302
303	/** \brief Default constructor.
304	*
305	* Initialises this->next=NULL and this->facetPtr=NULL
306	*/
307	gcone::gcone()
308	{
309	this->next=NULL;
310	this->prev=NULL;
311	this->facetPtr=NULL; //maybe this->facetPtr = new facet();
312	this->baseRing=currRing;
313	this->counter++;
314	this->UCN=this->counter;
315	this->numFacets=0;
316	this->ivIntPt=NULL;
317	}
318
319	/** \brief Constructor with ring and ideal
320	*
321	* This constructor takes the root ring and the root ideal as parameters and stores
322	* them in the private members gcone::rootRing and gcone::inputIdeal
323	* Since knowledge of the root ring is only needed when using reverse search,
324	* this constructor is not needed when using the "second" method
325	*/
326	gcone::gcone(ring r, ideal I)
327	{
328	this->next=NULL;
329	this->prev=NULL;
330	this->facetPtr=NULL;
331	this->rootRing=r;
332	this->inputIdeal=I;
333	this->baseRing=currRing;
334	this->counter++;
335	this->UCN=this->counter;
336	this->numFacets=0;
337	this->ivIntPt=NULL;
338	}
339
340	/** \brief Copy constructor
341	*
342	* Copies a cone, sets this->gcBasis to the flipped GB
343	* Call this only after a successful call to gcone::flip which sets facet::flipGB
344	*/
345	gcone::gcone(const gcone& gc, const facet &f)
346	{
347	this->next=NULL;
348	// this->prev=(gcone *)&gc; //comment in to get a tree
349	this->prev=NULL;
350	this->numVars=gc.numVars;
351	this->counter++;
352	this->UCN=this->counter;
353	this->facetPtr=NULL;
354	this->gcBasis=idCopy(f.flipGB);
355	this->inputIdeal=idCopy(this->gcBasis);
356	this->baseRing=rCopy(f.flipRing);
357	this->numFacets=0;
358	this->ivIntPt=NULL;
359	this->rootRing=NULL;
360	//rComplete(this->baseRing);
361	//rChangeCurrRing(this->baseRing);
362	}
363
364	/** \brief Default destructor
365	*/
366	gcone::~gcone()
367	{
368	if(this->gcBasis!=NULL)
369	idDelete((ideal *)&this->gcBasis);
370	if(this->inputIdeal!=NULL)
371	idDelete((ideal *)&this->inputIdeal);
372	if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
373	rDelete(this->rootRing);
374	//rDelete(this->baseRing);
375	facet *fAct;
376	facet *fDel;
377	/Delete the facet structure/
378	fAct=this->facetPtr;
379	fDel=fAct;
380	while(fAct!=NULL)
381	{
382	fDel=fAct;
383	fAct=fAct->next;
384	delete fDel;
385	}
386	//dd_FreeMatrix(this->ddFacets);
387	}
388
389	/** \brief Set the interior point of a cone */
390	void gcone::setIntPoint(intvec *iv)
391	{
392	this->ivIntPt=ivCopy(iv);
393	}
394
395	/** \brief Return the interior point */
396	intvec *gcone::getIntPoint()
397	{
398	return this->ivIntPt;
399	}
400
401	/** \brief Print the interior point */
402	void gcone::showIntPoint()
403	{
404	ivIntPt->show();
405	}
406
407	/** \brief Print facets
408	* This is mainly for debugging purposes. Usually called from within gdb
409	*/
410	void gcone::showFacets(short codim)
411	{
412	facet *f;
413	switch(codim)
414	{
415	case 1:
416	f = this->facetPtr;
417	break;
418	case 2:
419	f = this->facetPtr->codim2Ptr;
420	break;
421	}
422
423	intvec *iv;
424	while(f!=NULL)
425	{
426	iv = f->getFacetNormal();
427	cout << "(";
428	iv->show(1,0);
429	if(f->isFlippable==FALSE)
430	cout << ")* ";
431	else
432	cout << ") ";
433	f=f->next;
434	}
435	cout << endl;
436	}
437
438	/** For debugging purposes only */
439	void gcone::showSLA(facet &f)
440	{
441	facet *fAct;
442	fAct = &f;
443	facet *codim2Act;
444	codim2Act = fAct->codim2Ptr;
445	intvec *fNormal;
446	intvec *f2Normal;
447	cout << endl;
448	while(fAct!=NULL)
449	{
450	fNormal=fAct->getFacetNormal();
451	cout << "(";
452	fNormal->show(1,0);
453	if(fAct->isFlippable==TRUE)
454	cout << ") ";
455	else
456	cout << ")* ";
457	codim2Act = fAct->codim2Ptr;
458	cout << " Codim2: ";
459	while(codim2Act!=NULL)
460	{
461	f2Normal = codim2Act->getFacetNormal();
462	cout << "(";
463	f2Normal->show(1,0);
464	cout << ") ";
465	codim2Act = codim2Act->next;
466	}
467	cout << "UCN = " << fAct->getUCN() << endl;
468	fAct = fAct->next;
469	}
470	}
471
472	void gcone::idDebugPrint(ideal const &I)
473	{
474	int numElts=IDELEMS(I);
475	cout << "Ideal with " << numElts << " generators" << endl;
476	cout << "Leading terms: ";
477	for (int ii=0;ii<numElts;ii++)
478	{
479	pWrite0(pHead(I->m[ii]));
480	cout << ",";
481	}
482	cout << endl;
483	}
484
485	void gcone::invPrint(ideal const &I)
486	{
487	// int numElts=IDELEMS(I);
488	// cout << "inv = ";
489	// for(int ii=0;ii<numElts;ii++);
490	// {
491	// pWrite0(pHead(I->m[ii]));
492	// cout << ",";
493	// }
494	// cout << endl;
495	}
496
497	bool gcone::isMonomial(ideal const &I)
498	{
499	bool res = TRUE;
500	for(int ii=0;ii<IDELEMS(I);ii++)
501	{
502	if(pLength((poly)I->m[ii])>1)
503	{
504	res = FALSE;
505	break;
506	}
507	}
508	return res;
509	}
510
511	/** \brief Set gcone::numFacets */
512	void gcone::setNumFacets()
513	{
514	}
515
516	/** \brief Get gcone::numFacets */
517	int gcone::getNumFacets()
518	{
519	return this->numFacets;
520	}
521
522	int gcone::getUCN()
523	{
524	if( this!=NULL && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
525	return this->UCN;
526	else
527	return -1;
528	}
529	/** \brief Compute the normals of the cone
530	*
531	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
532	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
533	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
534	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
535	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
536	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
537	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
538	*
539	* Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
540	* an interior point of the cone.
541	*/
542	void gcone::getConeNormals(ideal const &I, bool compIntPoint)
543	{
544	#ifdef gfan_DEBUG
545	std::cout << "* Computing Inequalities... *" << std::endl;
546	#endif
547	//All variables go here - except ineq matrix and *v, see below
548	int lengthGB=IDELEMS(I); // # of polys in the groebner basis
549	int pCompCount; // # of terms in a poly
550	poly aktpoly;
551	int numvar = pVariables; // # of variables in a polynomial (or ring?)
552	int leadexp[numvar]; // dirty hack of exp.vects
553	int aktexp[numvar];
554	int cols,rows; // will contain the dimensions of the ineq matrix - deprecated by
555	dd_rowrange ddrows;
556	dd_colrange ddcols;
557	dd_rowset ddredrows; // # of redundant rows in ddineq
558	dd_rowset ddlinset; // the opposite
559	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
560	dd_NumberType ddnumb=dd_Integer; //Number type
561	dd_ErrorType dderr=dd_NoError; //
562	// End of var declaration
563	#ifdef gfan_DEBUG
564	// cout << "The Groebner basis has " << lengthGB << " elements" << endl;
565	// cout << "The current ring has " << numvar << " variables" << endl;
566	#endif
567	cols = numvar;
568
569	//Compute the # inequalities i.e. rows of the matrix
570	rows=0; //Initialization
571	for (int ii=0;ii<IDELEMS(I);ii++)
572	{
573	aktpoly=(poly)I->m[ii];
574	rows=rows+pLength(aktpoly)-1;
575	}
576	#ifdef gfan_DEBUG
577	// cout << "rows=" << rows << endl;
578	// cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
579	#endif
580	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
581	dd_set_global_constants();
582	ddrows=rows;
583	ddcols=cols;
584	dd_MatrixPtr ddineq; //Matrix to store the inequalities
585	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
586
587	// We loop through each g\in GB and compute the resulting inequalities
588	for (int i=0; i<IDELEMS(I); i++)
589	{
590	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
591	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
592	#ifdef gfan_DEBUG
593	// cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
594	#endif
595
596	int v=(int )omAlloc((numvar+1)*sizeof(int));
597	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
598
599	//Store leadexp for aktpoly
600	for (int kk=0;kk<numvar;kk++)
601	{
602	leadexp[kk]=v[kk+1];
603	//Since we need to know the difference of leadexp with the other expvects we do nothing here
604	//but compute the diff below
605	}
606
607
608	while (pNext(aktpoly)!=NULL) //move to next term until NULL
609	{
610	aktpoly=pNext(aktpoly);
611	pSetm(aktpoly); //doesn't seem to help anything
612	pGetExpV(aktpoly,v);
613
614	for (int kk=0;kk<numvar;kk++)
615	{
616	aktexp[kk]=v[kk+1];
617	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
618	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
619	}
620	aktmatrixrow=aktmatrixrow+1;
621	}//while
622	omFree(v);
623	} //for
624
625	#if false
626	/*Let's make the preprocessing here. This could already be done in the above for-loop,
627	* but for a start it is more convenient here.
628	* We check the necessary condition of FJT p.18
629	* Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
630	*/
631	ideal initialForm=idInit(IDELEMS(I),1);
632	intvec *gamma=new intvec(this->numVars);
633	int falseGammaCounter=0;
634	int *redRowsArray=NULL;
635	int num_alloc=0;
636	int num_elts=0;
637	for(int ii=0;ii<ddineq->rowsize;ii++)
638	{
639	double tmp;
640	for(int jj=1;jj<=this->numVars;jj++)
641	{
642	tmp=mpq_get_d(ddineq->matrix[ii][jj]);
643	(*gamma)[jj-1]=(int)tmp;
644	// (*gamma)[jj]=(ddineq->matrix[ii][jj]);
645	}
646	computeInv((ideal&)I,initialForm,*gamma);
647	//Create leading ideal
648	ideal L=idInit(IDELEMS(initialForm),1);
649	for(int jj=0;jj<IDELEMS(initialForm);jj++)
650	{
651	L->m[jj]=pCopy(pHead(initialForm->m[jj]));
652	}
653
654	LObject *P = new sLObject();
655	memset(P,0,sizeof(LObject));
656	// P->p1=L->m[0];
657	// P->p2=L->m[1];
658	// ksCreateSpoly(P);
659	// pWrite(P->p);
660	// cout << "gamma=";
661	// gamma->show(1,0);
662	// cout << endl;
663	// cout << "inv=";
664	// for(int qq=0;qq<IDELEMS(initialForm);qq++)
665	// {
666	// pWrite0(initialForm->m[qq]);
667	// cout << ",";
668	// }
669	// cout << endl;
670	for(int jj=0;jj<=IDELEMS(initialForm)-2;jj++)
671	{
672	bool isMaybeFacet=FALSE;
673	P->p1=initialForm->m[jj]; //build spolys of initialForm in_v
674	// cout << "P->p1=";
675	// pWrite0(P->p1);
676	// cout << endl;
677	for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
678	{
679	P->p2=initialForm->m[kk];
680	ksCreateSpoly(P);
681	/*cout << "P->p2=";
682	pWrite0(P->p2);
683	cout << endl;
684	cout << "spoly(p1,p2)=";
685	pWrite0(P->p);
686	cout << endl; */
687	if(P->p!=NULL) //spoly non zero=?
688	{
689	poly p=pInit();
690	poly q=pInit();
691	p=pCopy(P->p);
692	q=pHead(p); //Monomial q
693	isMaybeFacet=FALSE;
694	while(p!=NULL)
695	{
696	q=pHead(p);
697	// unsigned long sevSpoly=pGetShortExpVector(q);
698	// unsigned long not_sevL;
699	for(int ll=0;ll<IDELEMS(L);ll++)
700	{
701	// not_sevL=~pGetShortExpVector(L->m[ll]);//
702	//if(!(sevSpoly & not_sevL) && pLmDivisibleBy(L->m[ll],q) )//i.e. spoly is in L
703	if(pLmEqual(L->m[ll],q) \|\| pDivisibleBy(L->m[ll],q))
704	{
705	isMaybeFacet=TRUE;
706	break;
707	}
708	}
709	if(isMaybeFacet==TRUE)
710	{
711	// dd_MatrixRowRemove(&ddineq,ii);
712	// set_addelem(redRows,ii);
713	// dd_set_si(ddineq->matrix[ii][0],1);
714	// falseGammaCounter++;
715	// cout << "Removing gamma!" << endl;
716	break;//while(p!=NULL)
717	}
718	p=pNext(p);
719	}//while
720	if(isMaybeFacet==FALSE)
721	{
722	dd_set_si(ddineq->matrix[ii][0],1);
723	if(num_alloc==0)
724	num_alloc += 1;
725	else
726	num_alloc += 1;
727	void _tmp = realloc(redRowsArray,(num_allocsizeof(int)));
728	if(!_tmp)
729	WerrorS("Couldn't realloc memory\n");
730	redRowsArray = (int*)_tmp;
731	redRowsArray[num_elts]=ii;
732	num_elts++;
733	break;//for(int kk, since we have found one that is not in L
734	}
735	}//if(P->p!=NULL)
736	}
737	}//for
738	delete P;
739	//idDelete(L);
740	}//for(ii<ddineq-rowsize
741
742	int offset=0;//needed for correction of redRowsArray[ii]
743	for( int ii=0;ii<num_elts;ii++ )
744	{
745	dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
746	offset++;
747	}
748	free(redRowsArray);
749	#endif
750	// cout << "Found " << falseGammaCounter << " false in " << rows << endl;
751	//Maybe add another row to contain the constraints of the standard simplex?
752
753	#ifdef gfan_DEBUG
754	// cout << "The inequality matrix is" << endl;
755	// dd_WriteMatrix(stdout, ddineq);
756	#endif
757
758	// The inequalities are now stored in ddineq
759	// Next we check for superflous rows
760	// time_t canonicalizeTic, canonicalizeTac;
761	// time(&canonicalizeTic);
762	// ddredrows = dd_RedundantRows(ddineq, &dderr);
763	// if (dderr!=dd_NoError) // did an error occur?
764	// {
765	// dd_WriteErrorMessages(stderr,dderr); //if so tell us
766	// }
767	#ifdef gfan_DEBUG
768	// else
769	// {
770	// cout << "Redundant rows: ";
771	// set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
772	// }//if dd_Error
773	#endif
774	//Remove reduntant rows here!
775	//Necessary check here! C.f. FJT p18
776
777	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
778	// time(&canonicalizeTac);
779	// cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl;
780	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
781	ddcols = ddineq->colsize;
782
783	//ddCreateMatrix(ddrows,ddcols+1);
784	ddFacets = dd_CopyMatrix(ddineq);
785	#ifdef gfan_DEBUG
786	// cout << "Having removed redundancies, the normals now read:" << endl;
787	// dd_WriteMatrix(stdout,ddineq);
788	// cout << "Rows = " << ddrows << endl;
789	// cout << "Cols = " << ddcols << endl;
790	#endif
791
792	/Write the normals into class facet/
793	#ifdef gfan_DEBUG
794	// cout << "Creating list of normals" << endl;
795	#endif
796	/The pointer fRoot should be the return value of this function*/
797	//facet *fRoot = new facet(); //instantiate new facet
798	//fRoot->setUCN(this->getUCN());
799	//this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
800	facet *fAct; //pointer to active facet
801	//fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress!
802	//std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
803	int numNonFlip=0;
804	for (int kk = 0; kk<ddrows; kk++)
805	{
806	intvec *load = new intvec(this->numVars); //intvec to store a single facet normal that will then be stored via setFacetNormal
807	for (int jj = 1; jj <ddcols; jj++)
808	{
809	double foo;
810	foo = mpq_get_d(ddineq->matrix[kk][jj]);
811	#ifdef gfan_DEBUG
812	// std::cout << "fAct is " << foo << " at " << fAct << std::endl;
813	#endif
814	(*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
815	}//for (int jj = 1; jj <ddcols; jj++)
816
817	/Quick'n'dirty hack for flippability/
818	bool isFlip=FALSE;
819
820	for (int jj = 0; jj<load->length(); jj++)
821	{
822	intvec *ivCanonical = new intvec(load->length());
823	(*ivCanonical)[jj]=1;
824	// cout << "dotProd=" << dotProduct(load,ivCanonical) << endl;
825	if (dotProduct(load,ivCanonical)<0)
826	//if (ivMult(load,ivCanonical)<0)
827	{
828	isFlip=TRUE;
829	break; //URGHS
830	}
831	}
832	if (isFlip==FALSE)
833	{
834	this->numFacets++;
835	numNonFlip++;
836	if(this->numFacets==1)
837	{
838	facet *fRoot = new facet();
839	this->facetPtr = fRoot;
840	fAct = fRoot;
841
842	}
843	else
844	{
845	fAct->next = new facet();
846	fAct = fAct->next;
847	}
848	fAct->isFlippable=FALSE;
849	fAct->setFacetNormal(load);
850	fAct->setUCN(this->getUCN());
851	#ifdef gfan_DEBUG
852	cout << "Marking facet (";
853	load->show(1,0);
854	cout << ") as non flippable" << endl;
855	#endif
856	}
857	else
858	{ /Now load should be full and we can call setFacetNormal/
859	this->numFacets++;
860	if(this->numFacets==1)
861	{
862	facet *fRoot = new facet();
863	this->facetPtr = fRoot;
864	fAct = fRoot;
865	}
866	else
867	{
868	fAct->next = new facet();
869	fAct = fAct->next;
870	}
871	fAct->isFlippable=TRUE;
872	fAct->setFacetNormal(load);
873	fAct->setUCN(this->getUCN());
874	}//if (isFlippable==FALSE)
875	//delete load;
876	}//for (int kk = 0; kk<ddrows; kk++)
877
878	//In cases like I=<x-1,y-1> there are only non-flippable facets...
879	if(numNonFlip==this->numFacets)
880	{
881	cerr << "Only non-flippable facets. Terminating..." << endl;
882	}
883
884	/*
885	Now we should have a linked list containing the facet normals of those facets that are
886	-irredundant
887	-flipable
888	Adressing is done via *facetPtr
889	*/
890
891	if (compIntPoint==TRUE)
892	{
893	intvec *iv = new intvec(this->numVars);
894	dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
895	int jj=1;
896	for (int ii=0;ii<=this->numVars;ii++)
897	{
898	dd_set_si(posRestr->matrix[ii][jj],1);
899	jj++;
900	}
901	dd_MatrixAppendTo(&ddineq,posRestr);
902	interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets
903	this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt
904	//delete iv;
905	}
906
907	//Compute the number of facets
908	// wrong because ddineq->rowsize contains only irredundant facets. There can still be
909	// non-flippable ones!
910	//this->numFacets=ddineq->rowsize;
911
912	//Clean up but don't delete the return value! (Whatever it will turn out to be)
913	dd_FreeMatrix(ddineq);
914	set_free(ddredrows);
915	//free(ddnewpos);
916	//set_free(ddlinset);
917	//NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
918	//THIS SUCKS
919	//dd_free_global_constants();
920
921	}//method getConeNormals(ideal I)
922
923	/** \brief Compute the (codim-2)-facets of a given cone
924	* This method is used during noRevS
925	* Additionally we check whether the codim2-facet normal is strictly positive. Otherwise
926	* the facet is marked as non-flippable.
927	*/
928	void gcone::getCodim2Normals(gcone const &gc)
929	{
930	//this->facetPtr->codim2Ptr = new facet(2); //instantiate a (codim-2)-facet
931	facet *fAct;
932	fAct = this->facetPtr;
933	facet *codim2Act;
934	//codim2Act = this->facetPtr->codim2Ptr;
935
936	dd_set_global_constants();
937
938	dd_MatrixPtr ddineq,P,ddakt;
939	dd_rowset impl_linset, redset;
940	dd_ErrorType err;
941	dd_rowindex newpos;
942
943	//ddineq = facets2Matrix(gc); //get a matrix representation of the cone
944	ddineq = dd_CopyMatrix(gc.ddFacets);
945
946	/Now set appropriate linearity/
947	dd_PolyhedraPtr ddpolyh;
948	for (int ii=0; ii<this->numFacets; ii++)
949	{
950	ddakt = dd_CopyMatrix(ddineq);
951	ddakt->representation=dd_Inequality;
952	set_addelem(ddakt->linset,ii+1);
953	dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
954
955	#ifdef gfan_DEBUG
956	// cout << "Codim2 matrix"<<endl;
957	// dd_WriteMatrix(stdout,ddakt);
958	#endif
959	ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
960	P=dd_CopyGenerators(ddpolyh);
961	#ifdef gfan_DEBUG
962	// cout << "Codim2 facet:" << endl;
963	// dd_WriteMatrix(stdout,P);
964	// cout << endl;
965	#endif
966
967	/* We loop through each row of P
968	* normalize it by making all entries integer ones
969	* and add the resulting vector to the int matrix facet::codim2Facets
970	*/
971	for (int jj=1;jj<=P->rowsize;jj++)
972	{
973	fAct->numCodim2Facets++;
974	if(fAct->numCodim2Facets==1)
975	{
976	fAct->codim2Ptr = new facet(2);
977	codim2Act = fAct->codim2Ptr;
978	}
979	else
980	{
981	codim2Act->next = new facet(2);
982	codim2Act = codim2Act->next;
983	}
984	intvec *n = new intvec(this->numVars);
985	makeInt(P,jj,*n);
986	codim2Act->setFacetNormal(n);
987	delete n;
988	}
989	/*We check whether the facet spanned by the codim-2 facets
990	* intersects with the positive orthant. Otherwise we define this
991	* facet to be non-flippable
992	*/
993	intvec *iv_intPoint = new intvec(this->numVars);
994	dd_MatrixPtr shiftMatrix;
995	dd_MatrixPtr intPointMatrix;
996	shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
997	for(int kk=0;kk<this->numVars;kk++)
998	{
999	dd_set_si(shiftMatrix->matrix[kk][0],1);
1000	dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
1001	}
1002	intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
1003	//dd_WriteMatrix(stdout,intPointMatrix);
1004	interiorPoint(intPointMatrix,*iv_intPoint);
1005	for(int ll=0;ll<this->numVars;ll++)
1006	{
1007	if( (*iv_intPoint)[ll] < 0 )
1008	{
1009	fAct->isFlippable=FALSE;
1010	break;
1011	}
1012	}
1013	/End of check/
1014	fAct = fAct->next;
1015	dd_FreeMatrix(ddakt);
1016	// dd_FreeMatrix(ddineq);
1017	dd_FreePolyhedra(ddpolyh);
1018	delete iv_intPoint;
1019	}//while
1020	}
1021
1022	/** \brief Compute the Groebner Basis on the other side of a shared facet
1023	*
1024	* Implements algorithm 4.3.2 from Anders' thesis.
1025	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
1026	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
1027	* is parallel to \f$ leadexp(g) \f$
1028	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
1029	* Other possibilities include computing the rank of the matrix consisting of the vectors in question and
1030	* computing an interior point of the facet and taking all terms having the same weight with respect
1031	* to this interior point.
1032	*\param ideal, facet
1033	* Input: a marked,reduced Groebner basis and a facet
1034	*/
1035	void gcone::flip(ideal gb, facet *f) //Compute "the other side"
1036	{
1037	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
1038	fNormal = f->getFacetNormal(); //read this->fNormal;
1039	//#ifdef gfan_DEBUG
1040	//std::cout << "===" << std::endl;
1041	std::cout << "running gcone::flip" << std::endl;
1042	std::cout << "flipping UCN " << this->getUCN() << endl;
1043	// for(int ii=0;ii<IDELEMS(gb);ii++) //not very handy with large examples
1044	// {
1045	// pWrite((poly)gb->m[ii]);
1046	// }
1047	cout << "over facet (";
1048	fNormal->show(1,0);
1049	cout << ") with UCN " << f->getUCN();
1050	std::cout << std::endl;
1051	//#endif
1052	/1st step: Compute the initial ideal/
1053	/poly initialFormElement[IDELEMS(gb)];/ //array of #polys in GB to store initial form
1054	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
1055	poly aktpoly;
1056	/intvec check = new intvec(this->numVars);*/ //array to store the difference of LE and v
1057	computeInv(gb,initialForm,*fNormal);
1058	//NOTE The code below went into gcone::computeInv
1059	// for (int ii=0;ii<IDELEMS(gb);ii++)
1060	// {
1061	// aktpoly = (poly)gb->m[ii];
1062	// int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1063	// int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
1064	// pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
1065	// initialFormElement[ii]=pHead(aktpoly);
1066	//
1067	// while(pNext(aktpoly)!=NULL) //loop trough terms and check for parallelity/
1068	// {
1069	// aktpoly=pNext(aktpoly); //next term
1070	// pSetm(aktpoly);
1071	// pGetExpV(aktpoly,v);
1072	// /* Convert (int)v into (intvec)check */
1073	// for (int jj=0;jj<this->numVars;jj++)
1074	// {
1075	// //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
1076	// //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
1077	// (*check)[jj]=v[jj+1]-leadExpV[jj+1];
1078	// }
1079	// #ifdef gfan_DEBUG
1080	// // cout << "check=";
1081	// // check->show();
1082	// // cout << endl;
1083	// #endif
1084	// if (isParallel(check,fNormal)) //pass *check when
1085	// {
1086	// // cout << "Parallel vector found, adding to initialFormElement" << endl;
1087	// initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
1088	// }
1089	// }//while
1090	// #ifdef gfan_DEBUG
1091	// // cout << "Initial Form=";
1092	// // pWrite(initialFormElement[ii]);
1093	// // cout << "---" << endl;
1094	// #endif
1095	// //Now initialFormElement must be added to (ideal)initialForm /
1096	// initialForm->m[ii]=initialFormElement[ii];
1097	// omFree(leadExpV);
1098	// omFree(v);
1099	// }//for
1100	#ifdef gfan_DEBUG
1101	/* cout << "Initial ideal is: " << endl;
1102	idShow(initialForm);
1103	//f->printFlipGB();*/
1104	// cout << "===" << endl;
1105	#endif
1106	//delete check;
1107
1108	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
1109	/*Substep 2.1
1110	compute $G_{-\alpha}(in_v(I))
1111	see journal p. 66
1112	NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
1113	srcRing already has a weighted ordering
1114	*/
1115	ring srcRing=currRing;
1116	ring tmpRing;
1117
1118	if( (srcRing->order[0]!=ringorder_a))
1119	{
1120	tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
1121	}
1122	else
1123	{
1124	tmpRing=rCopy0(srcRing);
1125	int length=fNormal->length();
1126	int A=(int )omAlloc0(length*sizeof(int));
1127	for(int jj=0;jj<length;jj++)
1128	{
1129	A[jj]=-(*fNormal)[jj];
1130	}
1131	omFree(tmpRing->wvhdl[0]);
1132	tmpRing->wvhdl[0]=(int*)A;
1133	tmpRing->block1[0]=length;
1134	rComplete(tmpRing);
1135	}
1136	rChangeCurrRing(tmpRing);
1137
1138	//rWrite(currRing); cout << endl;
1139
1140	ideal ina;
1141	ina=idrCopyR(initialForm,srcRing);
1142	#ifndef NDEBUG
1143	#ifdef gfan_DEBUG
1144	cout << "ina=";
1145	idShow(ina); cout << endl;
1146	#endif
1147	#endif
1148	ideal H;
1149	//H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
1150	H=kStd(ina,NULL,testHomog,NULL); //This is \mathcal(G)_{>_-\alpha}(in_v(I))
1151	idSkipZeroes(H);
1152	#ifndef NDEBUG
1153	#ifdef gfan_DEBUG
1154	cout << "H="; idShow(H); cout << endl;
1155	#endif
1156	#endif
1157	/*Substep 2.2
1158	do the lifting and mark according to H
1159	*/
1160	rChangeCurrRing(srcRing);
1161	ideal srcRing_H;
1162	ideal srcRing_HH;
1163	srcRing_H=idrCopyR(H,tmpRing);
1164	#ifndef NDEBUG
1165	#ifdef gfan_DEBUG
1166	cout << "srcRing_H = ";
1167	idShow(srcRing_H); cout << endl;
1168	#endif
1169	#endif
1170	srcRing_HH=ffG(srcRing_H,this->gcBasis);
1171	#ifndef NDEBUG
1172	#ifdef gfan_DEBUG
1173	cout << "srcRing_HH = ";
1174	idShow(srcRing_HH); cout << endl;
1175	#endif
1176	#endif
1177	/*Substep 2.2.1
1178	* Mark according to G_-\alpha
1179	* Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
1180	* we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
1181	* represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
1182	* Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
1183	* compute the difference accordingly
1184	*/
1185	dd_set_global_constants();
1186	bool markingsAreCorrect=FALSE;
1187	dd_MatrixPtr intPointMatrix;
1188	int iPMatrixRows=0;
1189	dd_rowrange aktrow=0;
1190	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1191	{
1192	poly aktpoly=(poly)srcRing_HH->m[ii];
1193	iPMatrixRows = iPMatrixRows+pLength(aktpoly);
1194	}
1195	/* additionally one row for the standard-simplex and another for a row that becomes 0 during
1196	* construction of the differences
1197	*/
1198	intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1); //iPMatrixRows+10;
1199	intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
1200
1201	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1202	{
1203	markingsAreCorrect=FALSE; //crucial to initialise here
1204	poly aktpoly=srcRing_HH->m[ii];
1205	/Comparison of leading monomials is done via exponent vectors/
1206	for (int jj=0;jj<IDELEMS(H);jj++)
1207	{
1208	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
1209	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
1210	pGetExpV(aktpoly,src_ExpV);
1211	rChangeCurrRing(tmpRing); //this ring change is crucial!
1212	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
1213	rChangeCurrRing(srcRing);
1214	bool expVAreEqual=TRUE;
1215	for (int kk=1;kk<=this->numVars;kk++)
1216	{
1217	#ifdef gfan_DEBUG
1218	// cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
1219	#endif
1220	if (src_ExpV[kk]!=dst_ExpV[kk])
1221	{
1222	expVAreEqual=FALSE;
1223	}
1224	}
1225	//if (src_ExpV == dst_ExpV)
1226	if (expVAreEqual==TRUE)
1227	{
1228	markingsAreCorrect=TRUE; //everything is fine
1229	#ifdef gfan_DEBUG
1230	// cout << "correct markings" << endl;
1231	#endif
1232	}//if (pHead(aktpoly)==pHead(H->m[jj])
1233	// delete src_ExpV;
1234	// delete dst_ExpV;
1235	omFree(src_ExpV);
1236	omFree(dst_ExpV);
1237	}//for (int jj=0;jj<IDELEMS(H);jj++)
1238
1239	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1240	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
1241	if (markingsAreCorrect==TRUE)
1242	{
1243	pGetExpV(aktpoly,leadExpV);
1244	}
1245	else
1246	{
1247	rChangeCurrRing(tmpRing);
1248	pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
1249	rChangeCurrRing(srcRing);
1250	}
1251	/compute differences of the expvects/
1252	while (pNext(aktpoly)!=NULL)
1253	{
1254	/*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
1255	is not omitted when computing the differences*/
1256	if(markingsAreCorrect==TRUE)
1257	{
1258	aktpoly=pNext(aktpoly);
1259	pGetExpV(aktpoly,v);
1260	}
1261	else
1262	{
1263	pGetExpV(pHead(aktpoly),v);
1264	markingsAreCorrect=TRUE;
1265	}
1266
1267	for (int jj=0;jj<this->numVars;jj++)
1268	{
1269	/Store into ddMatrix/
1270	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
1271	}
1272	aktrow +=1;
1273	}
1274	// delete v;
1275	// delete leadExpV;
1276	omFree(v);
1277	omFree(leadExpV);
1278	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1279	/Now we add the constraint for the standard simplex/
1280	/NOTE:Might actually work without the standard simplex/
1281	dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
1282	for (int jj=1;jj<=this->numVars;jj++)
1283	{
1284	dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
1285	}
1286	//Let's make sure we compute interior points from the positive orthant
1287	dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1288	int jj=1;
1289	for (int ii=0;ii<this->numVars;ii++)
1290	{
1291	dd_set_si(posRestr->matrix[ii][jj],1);
1292	jj++;
1293	}
1294	dd_MatrixAppendTo(&intPointMatrix,posRestr);
1295	dd_FreeMatrix(posRestr);
1296	#ifdef gfan_DEBUG
1297	// dd_WriteMatrix(stdout,intPointMatrix);
1298	#endif
1299	intvec *iv_weight = new intvec(this->numVars);
1300	interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point
1301	dd_FreeMatrix(intPointMatrix);
1302	dd_free_global_constants();
1303
1304	/*Step 3
1305	turn the minimal basis into a reduced one
1306	*/
1307	//ring dstRing=rCopyAndAddWeight(tmpRing,iv_weight);
1308
1309	// NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
1310	// Thus:
1311	//ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
1312	//cout << "PLING" << endl;
1313	/*ring dstRing=rCopy0(srcRing);
1314	for (int ii=0;ii<this->numVars;ii++)
1315	{
1316	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
1317	}*/
1318	ring dstRing=rCopy0(tmpRing);
1319	int length=iv_weight->length();
1320	int A=(int )omAlloc0(length*sizeof(int));
1321	for(int jj=0;jj<length;jj++)
1322	{
1323	A[jj]=(*iv_weight)[jj];
1324	}
1325	dstRing->wvhdl[0]=(int*)A;
1326	rComplete(dstRing);
1327	rChangeCurrRing(dstRing);
1328	//rDelete(tmpRing);
1329	delete iv_weight;
1330	//#ifdef gfan_DEBUG
1331	rWrite(dstRing); cout << endl;
1332	//#endif
1333	ideal dstRing_I;
1334	dstRing_I=idrCopyR(srcRing_HH,srcRing);
1335	//dstRing_I=idrCopyR(inputIdeal,srcRing);
1336	//validOpts<1>=TRUE;
1337	#ifdef gfan_DEBUG
1338	//idShow(dstRing_I);
1339	#endif
1340	BITSET save=test;
1341	test\|=Sy_bit(OPT_REDSB);
1342	test\|=Sy_bit(OPT_REDTAIL);
1343	#ifdef gfan_DEBUG
1344	test\|=Sy_bit(6); //OPT_DEBUG
1345	#endif
1346	ideal tmpI;
1347	//NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
1348	//tmpI = idrCopyR(this->inputIdeal,this->baseRing);
1349	tmpI = dstRing_I;
1350	dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
1351	idNorm(dstRing_I);
1352	//kInterRed(dstRing_I);
1353	idSkipZeroes(dstRing_I);
1354	test=save;
1355	/End of step 3 - reduction/
1356
1357	f->setFlipGB(dstRing_I);//store the flipped GB
1358	f->flipRing=rCopy(dstRing); //store the ring on the other side
1359	//#ifdef gfan_DEBUG
1360	cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
1361	// f->printFlipGB();
1362	this->idDebugPrint(dstRing_I);
1363	cout << endl;
1364	//#endif
1365	rChangeCurrRing(srcRing); //return to the ring we started the computation of flipGB in
1366	// rDelete(dstRing);
1367	}//void flip(ideal gb, facet *f)
1368
1369	/** \brief Compute initial ideal
1370	* Compute the initial ideal in_v(G) wrt a (possible) facet normal
1371	* used in gcone::getFacetNormal in order to preprocess possible facet normals
1372	* and in gcone::flip for obvious reasons.
1373	*/
1374	void gcone::computeInv(ideal &gb, ideal &initialForm, intvec &fNormal)
1375	{
1376	intvec *check = new intvec(this->numVars);
1377	poly initialFormElement[IDELEMS(gb)];
1378	poly aktpoly;
1379
1380	for (int ii=0;ii<IDELEMS(gb);ii++)
1381	{
1382	aktpoly = (poly)gb->m[ii];
1383	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1384	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
1385	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
1386	initialFormElement[ii]=pHead(aktpoly);
1387
1388	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
1389	{
1390	aktpoly=pNext(aktpoly); //next term
1391	pSetm(aktpoly);
1392	pGetExpV(aktpoly,v);
1393	/* Convert (int)v into (intvec)check */
1394	for (int jj=0;jj<this->numVars;jj++)
1395	{
1396	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
1397	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
1398	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
1399	}
1400	#ifdef gfan_DEBUG
1401	// cout << "check=";
1402	// check->show();
1403	// cout << endl;
1404	#endif
1405	if (isParallel(check,fNormal)) //pass check when
1406	{
1407	// cout << "Parallel vector found, adding to initialFormElement" << endl;
1408	initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
1409	}
1410	}//while
1411	#ifdef gfan_DEBUG
1412	// cout << "Initial Form=";
1413	// pWrite(initialFormElement[ii]);
1414	// cout << "---" << endl;
1415	#endif
1416	/Now initialFormElement must be added to (ideal)initialForm /
1417	initialForm->m[ii]=initialFormElement[ii];
1418	omFree(leadExpV);
1419	omFree(v);
1420	}//for
1421	delete check;
1422	}
1423
1424	/** \brief Compute the remainder of a polynomial by a given ideal
1425	*
1426	* Compute \f$ f^{\mathcal{G}} \f$
1427	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
1428	* However, since we are only interested in the remainder, there is no need to
1429	* compute the factors \f$ a_i \f$
1430	*/
1431	//NOTE: Should be replaced by kNF or kNF2
1432	poly gcone::restOfDiv(poly const &f, ideal const &I)
1433	{
1434	// cout << "Entering restOfDiv" << endl;
1435	poly p=f;
1436	//pWrite(p);
1437	poly r=NULL; //The zero polynomial
1438	int ii;
1439	bool divOccured;
1440
1441	while (p!=NULL)
1442	{
1443	ii=1;
1444	divOccured=FALSE;
1445
1446	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
1447	{
1448	if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
1449	{
1450	poly step1,step2,step3;
1451	//cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
1452	step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
1453	//cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
1454	step2 = ppMult_qq(step1, I->m[ii-1]);
1455	step3 = pSub(pCopy(p), step2);
1456	//p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
1457	//pSetm(p);
1458	pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
1459	p=step3;
1460	//pWrite(p);
1461	divOccured=TRUE;
1462	}
1463	else
1464	{
1465	ii += 1;
1466	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
1467	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
1468	if (divOccured==FALSE)
1469	{
1470	//cout << "TICK 5" << endl;
1471	r=pAdd(pCopy(r),pHead(p));
1472	pSetm(r);
1473	pSort(r);
1474	//cout << "r="; pWrite(r); cout << endl;
1475
1476	if (pLength(p)!=1)
1477	{
1478	p=pSub(pCopy(p),pHead(p)); //Here it may occur that p=0 instead of p=NULL
1479	}
1480	else
1481	{
1482	p=NULL; //Hack to correct this situation
1483	}
1484	//cout << "p="; pWrite(p);
1485	}//if (divOccured==FALSE)
1486	}//while (p!=0)
1487	return r;
1488	}//poly restOfDiv(poly const &f, ideal const &I)
1489
1490	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
1491	*/
1492	//NOTE: use kNF or kNF2 instead of restOfDivision
1493	ideal gcone::ffG(ideal const &H, ideal const &G)
1494	{
1495	// cout << "Entering ffG" << endl;
1496	int size=IDELEMS(H);
1497	ideal res=idInit(size,1);
1498	poly temp1, temp2, temp3; //polys to temporarily store values for pSub
1499	for (int ii=0;ii<size;ii++)
1500	{
1501	// res->m[ii]=restOfDiv(H->m[ii],G);
1502	res->m[ii]=kNF(G, NULL,H->m[ii],0,0);
1503	temp1=H->m[ii];
1504	temp2=res->m[ii];
1505	temp3=pSub(temp1, temp2);
1506	res->m[ii]=temp3;
1507	//res->m[ii]=pSub(temp1,temp2); //buggy
1508	//pSort(res->m[ii]);
1509	//pSetm(res->m[ii]);
1510	//cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
1511	}
1512	return res;
1513	}
1514
1515	/** \brief Compute a Groebner Basis
1516	*
1517	* Computes the Groebner basis and stores the result in gcone::gcBasis
1518	*\param ideal
1519	*\return void
1520	*/
1521	void gcone::getGB(ideal const &inputIdeal)
1522	{
1523	BITSET save=test;
1524	test\|=Sy_bit(OPT_REDSB);
1525	test\|=Sy_bit(OPT_REDTAIL);
1526	ideal gb;
1527	gb=kStd(inputIdeal,NULL,testHomog,NULL);
1528	idNorm(gb);
1529	idSkipZeroes(gb);
1530	this->gcBasis=gb; //write the GB into gcBasis
1531	test=save;
1532	}//void getGB
1533
1534	/** \brief Compute the negative of a given intvec
1535	*/
1536	intvec gcone::ivNeg(const intvec iv)
1537	{
1538	intvec *res = new intvec(iv->length());
1539	res=ivCopy(iv);
1540	res = (int)-1;
1541	return res;
1542	}
1543
1544
1545	/** \brief Compute the dot product of two intvecs
1546	*
1547	*/
1548	int gcone::dotProduct(intvec const &iva, intvec const &ivb)
1549	{
1550	int res=0;
1551	for (int i=0;i<this->numVars;i++)
1552	{
1553	res = res+(iva[i]*ivb[i]);
1554	}
1555	return res;
1556	}//int dotProduct
1557
1558	/** \brief Check whether two intvecs are parallel
1559	*
1560	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1561	*/
1562	bool gcone::isParallel(intvec const &a, intvec const &b)
1563	{
1564	int lhs,rhs;
1565	bool res;
1566	lhs=dotProduct(a,b)*dotProduct(a,b);
1567	rhs=dotProduct(a,a)*dotProduct(b,b);
1568	//cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1569	if (lhs==rhs)
1570	{
1571	res = TRUE;
1572	}
1573	else
1574	{
1575	res = FALSE;
1576	}
1577	return res;
1578	}//bool isParallel
1579
1580	/** \brief Compute an interior point of a given cone
1581	* Result will be written into intvec iv.
1582	* Any rational point is automatically converted into an integer.
1583	*/
1584	void gcone::interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
1585	{
1586	dd_LPPtr lp,lpInt;
1587	dd_ErrorType err=dd_NoError;
1588	dd_LPSolverType solver=dd_DualSimplex;
1589	dd_LPSolutionPtr lpSol=NULL;
1590	dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior
1591	dd_rowindex ddnewpos;
1592	dd_NumberType numb;
1593	//M->representation=dd_Inequality;
1594	//M->objective-dd_LPMin; //Not sure whether this is needed
1595
1596	//NOTE: Make this n-dimensional!
1597	//dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1598
1599	/*NOTE: Leave the following line commented out!
1600	* Otherwise it will cause interior points that are not strictly positive on some examples
1601	*
1602	*/
1603	//dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1604	//if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1605	//cout << "Tick 2" << endl;
1606	//dd_WriteMatrix(stdout,M);
1607
1608	lp=dd_Matrix2LP(M, &err);
1609	if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}
1610	if (lp==NULL){cout << "LP is NULL" << endl;}
1611	#ifdef gfan_DEBUG
1612	// dd_WriteLP(stdout,lp);
1613	#endif
1614
1615	lpInt=dd_MakeLPforInteriorFinding(lp);
1616	if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1617	#ifdef gfan_DEBUG
1618	// dd_WriteLP(stdout,lpInt);
1619	#endif
1620
1621	dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1622	if (err!=dd_NoError)
1623	{
1624	cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1625	dd_WriteErrorMessages(stdout, err);
1626	}
1627
1628	//dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior
1629	if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1630
1631	//lpSol=dd_CopyLPSolution(lpInt);
1632	if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}
1633	#ifdef gfan_DEBUG
1634	cout << "Interior point: ";
1635	for (int ii=1; ii<(lpSol->d)-1;ii++)
1636	{
1637	dd_WriteNumber(stdout,lpSol->sol[ii]);
1638	}
1639	cout << endl;
1640	#endif
1641
1642	//NOTE The following strongly resembles parts of makeInt.
1643	//Maybe merge sometimes
1644	mpz_t kgV; mpz_init(kgV);
1645	mpz_set_str(kgV,"1",10);
1646	mpz_t den; mpz_init(den);
1647	mpz_t tmp; mpz_init(tmp);
1648	mpq_get_den(tmp,lpSol->sol[1]);
1649	for (int ii=1;ii<(lpSol->d)-1;ii++)
1650	{
1651	mpq_get_den(den,lpSol->sol[ii+1]);
1652	mpz_lcm(kgV,tmp,den);
1653	mpz_set(tmp, kgV);
1654	}
1655	mpq_t qkgV;
1656	mpq_init(qkgV);
1657	mpq_set_z(qkgV,kgV);
1658	for (int ii=1;ii<(lpSol->d)-1;ii++)
1659	{
1660	mpq_t product;
1661	mpq_init(product);
1662	mpq_mul(product,qkgV,lpSol->sol[ii]);
1663	iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
1664	mpq_clear(product);
1665	}
1666	#ifdef gfan_DEBUG
1667	// iv.show();
1668	// cout << endl;
1669	#endif
1670	mpq_clear(qkgV);
1671	mpz_clear(tmp);
1672	mpz_clear(den);
1673	mpz_clear(kgV);
1674
1675	dd_FreeLPSolution(lpSol);
1676	dd_FreeLPData(lpInt);
1677	dd_FreeLPData(lp);
1678	set_free(ddlinset);
1679	set_free(ddredrows);
1680
1681	}//void interiorPoint(dd_MatrixPtr const &M)
1682
1683	/** \brief Copy a ring and add a weighted ordering in first place
1684	*
1685	*/
1686	ring gcone::rCopyAndAddWeight(ring const &r, intvec const *ivw)
1687	{
1688	ring res=rCopy0(r);
1689	int jj;
1690
1691	omFree(res->order);
1692	res->order =(int )omAlloc0(4sizeof(int));
1693	omFree(res->block0);
1694	res->block0=(int )omAlloc0(4sizeof(int));
1695	omFree(res->block1);
1696	res->block1=(int )omAlloc0(4sizeof(int));
1697	omfree(res->wvhdl);
1698	res->wvhdl =(int *)omAlloc0(4sizeof(int**));
1699
1700	res->order[0]=ringorder_a;
1701	res->block0[0]=1;
1702	res->block1[0]=res->N;;
1703	res->order[1]=ringorder_dp; //basically useless, since that should never be used
1704	res->block0[1]=1;
1705	res->block1[1]=res->N;;
1706	res->order[2]=ringorder_C;
1707
1708	int length=ivw->length();
1709	int A=(int )omAlloc0(length*sizeof(int));
1710	for (jj=0;jj<length;jj++)
1711	{
1712	A[jj]=(*ivw)[jj];
1713	}
1714	res->wvhdl[0]=(int *)A;
1715	res->block1[0]=length;
1716
1717	rComplete(res);
1718	return res;
1719	}//rCopyAndAdd
1720
1721	ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1722	{
1723	ring res=rCopy0(currRing);
1724	rComplete(res);
1725	rSetWeightVec(res,(int64*)ivw);
1726	//rChangeCurrRing(rnew);
1727	return res;
1728	}
1729
1730	/** \brief Checks whether a given facet is a search facet
1731	* Determines whether a given facet of a cone is the search facet of a neighbouring cone
1732	* This is done in the following way:
1733	* We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1734	* that is first crossed during the generic walk.
1735	* We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1736	* If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1737	*/
1738	bool gcone::isSearchFacet(gcone &gcTmp, facet *testfacet)
1739	{
1740	ring actRing=currRing;
1741	facet facetPtr=(facet)gcTmp.facetPtr;
1742	facet fMin=new facet(facetPtr); //Pointer to the "minimal" facet
1743	//facet *fMin = new facet(tmpcone.facetPtr);
1744	//fMin=tmpcone.facetPtr; //Initialise to first facet of tmpcone
1745	facet *fAct; //Ptr to alpha_i
1746	facet *fCmp; //Ptr to alpha_j
1747	fAct = fMin;
1748	fCmp = fMin->next;
1749
1750	rChangeCurrRing(this->rootRing); //because we compare the monomials in the rootring
1751	poly p=pInit();
1752	poly q=pInit();
1753	intvec *alpha_i = new intvec(this->numVars);
1754	intvec *alpha_j = new intvec(this->numVars);
1755	intvec *sigma = new intvec(this->numVars);
1756	sigma=gcTmp.getIntPoint();
1757
1758	int u=(int )omAlloc((this->numVars+1)*sizeof(int));
1759	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1760	u[0]=0; v[0]=0;
1761	int weight1,weight2;
1762	while(fAct->next->next!=NULL) //NOTE this is ugly. Can it be done without fCmp?
1763	{
1764	/* Get alpha_i and alpha_{i+1} */
1765	alpha_i=fAct->getFacetNormal();
1766	alpha_j=fCmp->getFacetNormal();
1767	#ifdef gfan_DEBUG
1768	alpha_i->show();
1769	alpha_j->show();
1770	#endif
1771	/Compute the dot product of sigma and alpha_{i,j}/
1772	weight1=dotProduct(sigma,alpha_i);
1773	weight2=dotProduct(sigma,alpha_j);
1774	#ifdef gfan_DEBUG
1775	cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1776	#endif
1777	/Adjust alpha_i and alpha_i+1 accordingly/
1778	for(int ii=1;ii<=this->numVars;ii++)
1779	{
1780	u[ii]=weight1(alpha_i)[ii-1];
1781	v[ii]=weight2(alpha_j)[ii-1];
1782	}
1783
1784	/Now p_weight and q_weight need to be compared as exponent vectors/
1785	pSetCoeff0(p,nInit(1));
1786	pSetCoeff0(q,nInit(1));
1787	pSetExpV(p,u);
1788	pSetm(p);
1789	pSetExpV(q,v);
1790	pSetm(q);
1791	#ifdef gfan_DEBUG
1792	pWrite(p);pWrite(q);
1793	#endif
1794	/*We want to check whether x^p < x^q
1795	=> want to check for return value 1 */
1796	if (pLmCmp(p,q)==1) //i.e. x^q is smaller
1797	{
1798	fMin=fCmp;
1799	fAct=fMin;
1800	fCmp=fCmp->next;
1801	}
1802	else
1803	{
1804	//fAct=fAct->next;
1805	if(fCmp->next!=NULL)
1806	{
1807	fCmp=fCmp->next;
1808	}
1809	else
1810	{
1811	fAct=fAct->next;
1812	}
1813	}
1814	//fAct=fAct->next;
1815	}//while(fAct.facetPtr->next!=NULL)
1816	delete alpha_i,alpha_j,sigma;
1817
1818	/If testfacet was minimal then fMin should still point there /
1819
1820	//if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))
1821	#ifdef gfan_DEBUG
1822	cout << "Checking for parallelity" << endl <<" fMin is";
1823	fMin->printNormal();
1824	cout << "testfacet is ";
1825	testfacet->printNormal();
1826	cout << endl;
1827	#endif
1828	if (fMin==gcTmp.facetPtr)
1829	//if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
1830	//if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
1831	{
1832	cout << "Parallel" << endl;
1833	rChangeCurrRing(actRing);
1834	//delete alpha_min, test;
1835	return TRUE;
1836	}
1837	else
1838	{
1839	cout << "Not parallel" << endl;
1840	rChangeCurrRing(actRing);
1841	//delete alpha_min, test;
1842	return FALSE;
1843	}
1844	}//bool isSearchFacet
1845
1846	/** \brief Check for equality of two intvecs
1847	*/
1848	bool gcone::areEqual(intvec const &a, intvec const &b)
1849	{
1850	bool res=TRUE;
1851	for(int ii=0;ii<this->numVars;ii++)
1852	{
1853	if(a[ii]!=b[ii])
1854	{
1855	res=FALSE;
1856	break;
1857	}
1858	}
1859	return res;
1860	}
1861
1862	/** \brief The reverse search algorithm
1863	*/
1864	void gcone::reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1865	{
1866	facet *fAct=new facet();
1867	fAct = gcAct->facetPtr;
1868
1869	while(fAct!=NULL)
1870	{
1871	cout << "======================"<< endl;
1872	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1873	gcone gcTmp = new gcone(gcAct);
1874	//idShow(gcTmp->gcBasis);
1875	gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1876	#ifdef gfan_DEBUG
1877	facet *f = new facet();
1878	f=gcTmp->facetPtr;
1879	while(f!=NULL)
1880	{
1881	f->printNormal();
1882	f=f->next;
1883	}
1884	#endif
1885	gcTmp->showIntPoint();
1886	/recursive part goes gere/
1887	if (isSearchFacet(gcTmp,(facet)gcAct->facetPtr))
1888	{
1889	gcAct->next=gcTmp;
1890	cout << "PING"<< endl;
1891	reverseSearch(gcTmp);
1892	}
1893	else
1894	{
1895	delete gcTmp;
1896	/NOTE remove fAct from linked list. It's no longer needed/
1897	}
1898	/recursion ends/
1899	fAct = fAct->next;
1900	}//while(fAct->next!=NULL)
1901	}//reverseSearch
1902
1903	/** \brief The new method of Markwig and Jensen
1904	* Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1905	* In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1906	* Keep a list of facets as a linked list containing an intvec and an integer matrix.
1907	* Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1908	* soon as we compute this facet again. Comparison of facets is done by...
1909	*/
1910	void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
1911	{
1912	facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
1913	facet *SearchListAct;
1914	SearchListAct = NULL;
1915	//SearchListAct = SearchListRoot;
1916
1917	gcone *gcAct;
1918	gcAct = &gcRoot;
1919	gcone *gcPtr; //Pointer to end of linked list of cones
1920	gcPtr = &gcRoot;
1921	gcone *gcNext; //Pointer to next cone to be visited
1922	gcNext = &gcRoot;
1923	gcone *gcHead;
1924	gcHead = &gcRoot;
1925
1926	facet *fAct;
1927	fAct = gcAct->facetPtr;
1928
1929	ring rAct;
1930	rAct = currRing;
1931
1932	int UCNcounter=gcAct->getUCN();
1933
1934	/* Use pure SLA or writeCone2File */
1935	// enum heuristic {
1936	// ram,
1937	// hdd
1938	// };
1939	// heuristic h;
1940	// h=hdd;
1941
1942	#ifdef gfan_DEBUG
1943	cout << "NoRevs" << endl;
1944	cout << "Facets are:" << endl;
1945	gcAct->showFacets();
1946	#endif
1947
1948	gcAct->getCodim2Normals(*gcAct);
1949
1950	//Compute unique representation of codim-2-facets
1951	gcAct->normalize();
1952
1953	/* Make a copy of the facet list of first cone
1954	Since the operations getCodim2Normals and normalize affect the facets
1955	we must not memcpy them before these ops!
1956	*/
1957	intvec *fNormal;// = new intvec(this->numVars);
1958	intvec *f2Normal;// = new intvec(this->numVars);
1959	facet *codim2Act; codim2Act = NULL;
1960	facet *sl2Root; //sl2Root = new facet(2);
1961	facet *sl2Act; //sl2Act = sl2Root;
1962
1963	for(int ii=0;ii<this->numFacets;ii++)
1964	{
1965	//only copy flippable facets! NOTE: We do not need flipRing or any such information.
1966	if(fAct->isFlippable==TRUE)
1967	{
1968	fNormal = fAct->getFacetNormal();
1969	if( ii==0 \|\| (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
1970	{
1971	SearchListAct = SearchListRoot;
1972	//memcpy(SearchListAct,fAct,sizeof(facet));
1973	}
1974	else
1975	{
1976	SearchListAct->next = new facet();
1977	SearchListAct = SearchListAct->next;
1978	//memcpy(SearchListAct,fAct,sizeof(facet));
1979	}
1980	SearchListAct->setFacetNormal(fNormal);
1981	SearchListAct->setUCN(this->getUCN());
1982	SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
1983	SearchListAct->isFlippable=TRUE;
1984	//Copy codim2-facets
1985	codim2Act=fAct->codim2Ptr;
1986	SearchListAct->codim2Ptr = new facet(2);
1987	sl2Root = SearchListAct->codim2Ptr;
1988	sl2Act = sl2Root;
1989	//while(codim2Act!=NULL)
1990	for(int jj=0;jj<fAct->numCodim2Facets;jj++)
1991	{
1992	f2Normal = codim2Act->getFacetNormal();
1993	if(jj==0)
1994	{
1995	sl2Act = sl2Root;
1996	sl2Act->setFacetNormal(f2Normal);
1997	}
1998	else
1999	{
2000	sl2Act->next = new facet(2);
2001	sl2Act = sl2Act->next;
2002	sl2Act->setFacetNormal(f2Normal);
2003	}
2004	codim2Act = codim2Act->next;
2005	}
2006	fAct = fAct->next;
2007	}//if(fAct->isFlippable==TRUE)
2008	else {fAct = fAct->next;}
2009	}//End of copying facets into SLA
2010
2011	SearchListAct = SearchListRoot; //Set to beginning of list
2012	/Make SearchList doubly linked/
2013	while(SearchListAct!=NULL)
2014	{
2015	if(SearchListAct->next!=NULL)
2016	{
2017	SearchListAct->next->prev = SearchListAct;
2018	}
2019	SearchListAct = SearchListAct->next;
2020	}
2021	SearchListAct = SearchListRoot; //Set to beginning of List
2022
2023	fAct = gcAct->facetPtr;
2024	//NOTE Disabled until code works as expected. Reenabled 2.11.2009
2025	gcAct->writeConeToFile(*gcAct);
2026
2027	/End of initialisation/
2028	fAct = SearchListAct;
2029	/*2nd step
2030	Choose a facet from fListPtr, flip it and forget the previous cone
2031	We always choose the first facet from fListPtr as facet to be flipped
2032	*/
2033	time_t tic, tac;
2034	while((SearchListAct!=NULL))// && counter<10)
2035	{//NOTE See to it that the cone is only changed after ALL facets have been flipped!
2036	fAct = SearchListAct;
2037
2038	while(fAct!=NULL)
2039	{ //Since SLA should only contain flippables there should be no need to check for that
2040
2041	gcAct->flip(gcAct->gcBasis,fAct);
2042
2043
2044	ring rTmp=rCopy(fAct->flipRing);
2045	rComplete(rTmp);
2046	rChangeCurrRing(rTmp);
2047	gcone gcTmp = new gcone::gcone(gcAct,*fAct);
2048	/* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
2049	* Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
2050	* NOTE: Can this cause trouble with the destructor?
2051	* Answer: Yes it bloody well does!
2052	* However I'd like to delete this data here, since if we have a cone with many many facets it
2053	* should pay to get rid of the flipGB as soon as possible.
2054	* Destructor must be equipped with necessary checks.
2055	*/
2056	// idDelete((ideal *)&fAct->flipGB);
2057	// rDelete(fAct->flipRing);
2058
2059	gcTmp->getConeNormals(gcTmp->gcBasis, FALSE);
2060	gcTmp->getCodim2Normals(*gcTmp);
2061	gcTmp->normalize();
2062	#ifdef gfan_DEBUG
2063	gcTmp->showFacets(1);
2064	#endif
2065	if(gfanHeuristic==1)
2066	{
2067	gcTmp->writeConeToFile(*gcTmp);
2068	}
2069	/add facets to SLA here/
2070	SearchListRoot=gcTmp->enqueueNewFacets(*SearchListRoot);
2071	#ifdef gfan_DEBUG
2072	if(SearchListRoot!=NULL)
2073	gcTmp->showSLA(*SearchListRoot);
2074	#endif
2075	rChangeCurrRing(gcAct->baseRing);
2076	//rDelete(rTmp);
2077	//doubly linked for easier removal
2078	gcTmp->prev = gcPtr;
2079	gcPtr->next=gcTmp;
2080	gcPtr=gcPtr->next;
2081	if(fAct->getUCN() == fAct->next->getUCN())
2082	{
2083	fAct=fAct->next;
2084	}
2085	else
2086	break;
2087	}//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
2088	//NOTE Now all neighbouring cones of gcAct have been computed, so we may delete gcAct
2089	gcone *deleteMarker;
2090	deleteMarker=gcAct;
2091	// delete deleteMarker;
2092	// deleteMarker=NULL;
2093	//Search for cone with smallest UCN
2094	gcNext = gcHead;
2095	while(gcNext!=NULL && SearchListRoot!=NULL && gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && gcNext!=(gcone * const)0xfbfbfbfb)
2096	{
2097	/*NOTE if gcNext->getUCN is smaller than SearchListRoot->getUCN it follows, that the cone
2098	gcNext will not be used in any more computation. So -> delete
2099	*/
2100	if (gcNext->getUCN() < SearchListRoot->getUCN() )
2101	{
2102	//idDelete((ideal *)&gcNext->gcBasis); //at least drop huge grÃ¶bner bases
2103	if(gcNext == gcHead)
2104	{
2105	deleteMarker = gcHead;
2106	gcHead = gcNext->next;
2107	//gcNext->next->prev = NULL;
2108	}
2109	else
2110	{
2111	deleteMarker = gcNext;
2112	gcNext->prev->next = gcNext->next;
2113	gcNext->next->prev = gcNext->prev;
2114	}
2115	// gcNext = gcNext->next;
2116	// delete deleteMarker;
2117	// deleteMarker = NULL;
2118	}
2119	else if( gcNext->getUCN() == SearchListRoot->getUCN() )
2120	{//NOTE: Implement heuristic to be used!
2121	if(gfanHeuristic==0)
2122	{
2123	gcAct = gcNext;
2124	rAct=rCopy(gcAct->baseRing);
2125	rComplete(rAct);
2126	rChangeCurrRing(rAct);
2127	break;
2128	}
2129	else if(gfanHeuristic==1)
2130	{
2131	//Read st00f from file
2132	gcAct = gcNext;
2133	//implant the GB into gcAct
2134	readConeFromFile(gcNext->getUCN(), gcAct);
2135	rAct=rCopy(gcAct->baseRing);
2136	rComplete(rAct);
2137	rChangeCurrRing(rAct);
2138	break;
2139	}
2140	}
2141	gcNext = gcNext->next;
2142	// if(deleteMarker!=NULL && deleteMarker!=(gcone *)0xfbfbfbfbfbfbfbfb)
2143	if(this->getUCN()!=1) //workaround to avoid unpredictable behaviour towards end of program
2144	delete deleteMarker;
2145	// deleteMarker = NULL;
2146	// gcNext = gcNext->next;
2147	}
2148	UCNcounter++;
2149	//SearchListAct = SearchListAct->next;
2150	//SearchListAct = fAct->next;
2151	SearchListAct = SearchListRoot;
2152	}
2153	cout << endl << "Found " << counter << " cones - terminating" << endl;
2154
2155	//NOTE Hm, comment in and get a crash for free...
2156	//dd_free_global_constants();
2157	//gc.writeConeToFile(gc);
2158
2159
2160	//fAct = fListPtr;
2161	//gcone *gcTmp = new gcone(gc); //copy
2162	//gcTmp->flip(gcTmp->gcBasis,fAct);
2163	//NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
2164
2165	}//void noRevS(gcone &gc)
2166
2167
2168	/** \brief Make a set of rational vectors into integers
2169	*
2170	* We compute the lcm of the denominators and multiply with this to get integer values.
2171	* \param dd_MatrixPtr,intvec
2172	*/
2173	void gcone::makeInt(dd_MatrixPtr const &M, int const line, intvec &n)
2174	{
2175	mpz_t denom[this->numVars];
2176	for(int ii=0;ii<this->numVars;ii++)
2177	{
2178	mpz_init_set_str(denom[ii],"0",10);
2179	}
2180
2181	mpz_t kgV,tmp;
2182	mpz_init(kgV);
2183	mpz_init(tmp);
2184
2185	for (int ii=0;ii<(M->colsize)-1;ii++)
2186	{
2187	mpz_t z;
2188	mpz_init(z);
2189	mpq_get_den(z,M->matrix[line-1][ii+1]);
2190	mpz_set( denom[ii], z);
2191	mpz_clear(z);
2192	}
2193
2194	/Compute lcm of the denominators/
2195	mpz_set(tmp,denom[0]);
2196	for (int ii=0;ii<(M->colsize)-1;ii++)
2197	{
2198	mpz_lcm(kgV,tmp,denom[ii]);
2199	mpz_set(tmp,kgV);
2200	}
2201
2202	/Multiply the nominators by kgV/
2203	mpq_t qkgV,res;
2204	mpq_init(qkgV);
2205	mpq_set_str(qkgV,"1",10);
2206	mpq_canonicalize(qkgV);
2207
2208	mpq_init(res);
2209	mpq_set_str(res,"1",10);
2210	mpq_canonicalize(res);
2211
2212	mpq_set_num(qkgV,kgV);
2213
2214	// mpq_canonicalize(qkgV);
2215	for (int ii=0;ii<(M->colsize)-1;ii++)
2216	{
2217	mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
2218	n[ii]=(int)mpz_get_d(mpq_numref(res));
2219	}
2220	//mpz_clear(denom[this->numVars]);
2221	mpz_clear(kgV);
2222	mpq_clear(qkgV); mpq_clear(res);
2223
2224	}
2225	/**
2226	* For all codim-2-facets we compute the gcd of the components of the facet normal and
2227	* divide it out. Thus we get a normalized representation of each
2228	* (codim-2)-facet normal, i.e. a primitive vector
2229	*/
2230	void gcone::normalize()
2231	{
2232	int ggT=1;
2233	facet *fAct;
2234	facet *codim2Act;
2235	fAct = this->facetPtr;
2236	codim2Act = fAct->codim2Ptr;
2237	intvec *n = new intvec(this->numVars);
2238
2239	//while(codim2Act->next!=NULL)
2240	while(fAct!=NULL)
2241	{
2242	while(codim2Act!=NULL)
2243	{
2244	n=codim2Act->getFacetNormal();
2245	for(int ii=0;ii<this->numVars;ii++)
2246	{
2247	ggT = intgcd(ggT,(*n)[ii]);
2248	}
2249	for(int ii=0;ii<this->numVars;ii++)
2250	{
2251	(n)[ii] = ((n)[ii])/ggT;
2252	}
2253	codim2Act->setFacetNormal(n);
2254	codim2Act = codim2Act->next;
2255	}
2256	fAct = fAct->next;
2257	}
2258	delete n;
2259
2260	}
2261	/** \brief Enqueue new facets into the searchlist
2262	* The searchlist (SLA for short) is implemented as a FIFO
2263	* Checks are done as follows:
2264	* 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
2265	* 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
2266	* facet from SLA and do not add fAct.
2267	* It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
2268	* can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA.
2269	* Takes ptr to search list root
2270	* Returns a pointer to new first element of Searchlist
2271	*/
2272	//void enqueueNewFacets(facet &f)
2273	facet * gcone::enqueueNewFacets(facet &f)
2274	{
2275	facet *slHead;
2276	slHead = &f;
2277	facet *slAct; //called with f=SearchListRoot
2278	slAct = &f;
2279	facet *slEnd; //Pointer to end of SLA
2280	slEnd = &f;
2281	facet *slEndStatic; //marks the end before a new facet is added
2282	facet *fAct;
2283	fAct = this->facetPtr;
2284	facet *codim2Act;
2285	codim2Act = this->facetPtr->codim2Ptr;
2286	facet *sl2Act;
2287	sl2Act = f.codim2Ptr;
2288	/** Pointer to a facet that will be deleted */
2289	facet *deleteMarker;
2290	deleteMarker = NULL;
2291
2292	/** Flag to indicate a facet that should be added to SLA*/
2293	bool doNotAdd=FALSE;
2294	/** \brief Flag to mark a facet that might be added
2295	* The following case may occur:
2296	* A facet fAct is found to be parallel but not equal to the current facet slAct from SLA.
2297	* This does however not mean that there does not exist a facet behind the current slAct that is actually equal
2298	* to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
2299	* FALSE and the according slAct is deleted.
2300	* If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
2301	*/
2302	volatile bool maybe=FALSE;
2303	/**A facet was removed, lengthOfSearchlist-- thus we must not rely on
2304	* if(notParallelCtr==lengthOfSearchList) but rather
2305	* if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
2306	*/
2307	volatile bool removalOccured=FALSE;
2308	int ctr=0; //encountered equalities in SLA
2309	int notParallelCtr=0;
2310	int lengthOfSearchList=1;
2311	while(slEnd->next!=NULL)
2312	{
2313	slEnd=slEnd->next;
2314	lengthOfSearchList++;
2315	}
2316	slEndStatic = slEnd;
2317	/1st step: compare facetNormals/
2318	intvec *fNormal=NULL; //= new intvec(this->numVars);
2319	intvec *f2Normal=NULL; //= new intvec(this->numVars);
2320	intvec *slNormal=NULL; //= new intvec(this->numVars);
2321	intvec *sl2Normal=NULL; //= new intvec(this->numVars);
2322
2323	while(fAct!=NULL)
2324	{
2325	if(fAct->isFlippable==TRUE)
2326	{
2327	maybe=FALSE;
2328	doNotAdd=TRUE;
2329	fNormal=fAct->getFacetNormal();
2330	slAct = slHead;
2331	notParallelCtr=0;
2332	// delete deleteMarker;
2333	// deleteMarker=NULL;
2334	/*If slAct==NULL and fAct!=NULL
2335	we just copy all remaining facets into SLA*/
2336	if(slAct==NULL)
2337	{
2338	facet *fCopy;
2339	fCopy = fAct;
2340	while(fCopy!=NULL)
2341	{
2342	if(slAct==NULL)
2343	{
2344	slAct = new facet();
2345	slHead = slAct;
2346	}
2347	else
2348	{
2349	slAct->next = new facet();
2350	slAct = slAct->next;
2351	}
2352	slAct->setFacetNormal(fAct->getFacetNormal());
2353	slAct->setUCN(fAct->getUCN());
2354	slAct->isFlippable=fAct->isFlippable;
2355	facet *f2Copy;
2356	f2Copy = fCopy->codim2Ptr;
2357	sl2Act = slAct->codim2Ptr;
2358	while(f2Copy!=NULL)
2359	{
2360	if(sl2Act==NULL)
2361	{
2362	sl2Act = new facet(2);
2363	slAct->codim2Ptr = sl2Act;
2364	}
2365	else
2366	{
2367	facet *marker;
2368	marker = sl2Act;
2369	sl2Act->next = new facet(2);
2370	sl2Act = sl2Act->next;
2371	sl2Act->prev = marker;
2372	}
2373	sl2Act->setFacetNormal(f2Copy->getFacetNormal());
2374	sl2Act->setUCN(f2Copy->getUCN());
2375	f2Copy = f2Copy->next;
2376	}
2377	fCopy = fCopy->next;
2378	}
2379	break;
2380	}
2381	/End of dumping into SLA/
2382	while(slAct!=NULL)
2383	//while(slAct!=slEndStatic->next)
2384	{
2385	// if(deleteMarker!=NULL)
2386	// {
2387	// delete deleteMarker;
2388	// deleteMarker=NULL;
2389	// }
2390	removalOccured=FALSE;
2391	slNormal = slAct->getFacetNormal();
2392	#ifdef gfan_DEBUG
2393	cout << "Checking facet (";
2394	fNormal->show(1,1);
2395	cout << ") against (";
2396	slNormal->show(1,1);
2397	cout << ")" << endl;
2398	#endif
2399	if(!isParallel(fNormal, slNormal))
2400	{
2401	notParallelCtr++;
2402	// slAct = slAct->next;
2403	}
2404	else //fN and slN are parallel
2405	{
2406	//We check whether the codim-2-facets coincide
2407	codim2Act = fAct->codim2Ptr;
2408	ctr=0;
2409	while(codim2Act!=NULL)
2410	{
2411	f2Normal = codim2Act->getFacetNormal();
2412	//sl2Act = f.codim2Ptr;
2413	sl2Act = slAct->codim2Ptr;
2414	while(sl2Act!=NULL)
2415	{
2416	sl2Normal = sl2Act->getFacetNormal();
2417	if(areEqual(f2Normal, sl2Normal))
2418	ctr++;
2419	sl2Act = sl2Act->next;
2420	}//while(sl2Act!=NULL)
2421	codim2Act = codim2Act->next;
2422	}//while(codim2Act!=NULL)
2423	if(ctr==fAct->numCodim2Facets) //facets ARE equal
2424	{
2425	#ifdef gfan_DEBUG
2426	cout << "Removing facet (";
2427	slNormal->show(1,0);
2428	cout << ") from SLA:" << endl;
2429	fAct->fDebugPrint();
2430	slAct->fDebugPrint();
2431	#endif
2432	removalOccured=TRUE;
2433	slAct->isFlippable=FALSE;
2434	doNotAdd=TRUE;
2435	maybe=FALSE;
2436	if(slAct==slHead) //We want to delete the first element of SearchList
2437	{
2438	deleteMarker = slHead;
2439	slHead = slAct->next;
2440	if(slHead!=NULL)
2441	slHead->prev = NULL;
2442	//set a bool flag to mark slAct as to be deleted
2443	}//NOTE find a way to delete without affecting slAct = slAct->next
2444	else if(slAct==slEndStatic)
2445	{
2446	deleteMarker = slAct;
2447	if(slEndStatic->next==NULL)
2448	{
2449	slEndStatic = slEndStatic->prev;
2450	slEndStatic->next = NULL;
2451	slEnd = slEndStatic;
2452	}
2453	else //we already added a facet after slEndStatic
2454	{
2455	slEndStatic->prev->next = slEndStatic->next;
2456	slEndStatic->next->prev = slEndStatic->prev;
2457	slEndStatic = slEndStatic->prev;
2458	//slEnd = slEndStatic;
2459	}
2460	}
2461	else
2462	{
2463	deleteMarker = slAct;
2464	slAct->prev->next = slAct->next;
2465	slAct->next->prev = slAct->prev;
2466	}
2467
2468	//update lengthOfSearchList
2469	lengthOfSearchList--;
2470	//delete slAct;
2471	//slAct=NULL;
2472	//slAct = slAct->next; //not needed, since facets are equal
2473	//delete deleteMarker;
2474	//deleteMarker=NULL;
2475	//fAct = fAct->next;
2476	break;
2477	}//if(ctr==fAct->numCodim2Facets)
2478	else //facets are parallel but NOT equal. But this does not imply that there
2479	//is no other facet later in SLA that might be equal.
2480	{
2481	maybe=TRUE;
2482	// if(slAct->next==NULL && maybe==TRUE)
2483	// {
2484	// doNotAdd=FALSE;
2485	// slAct = slAct->next;
2486	// break;
2487	// }
2488	// else
2489	// slAct=slAct->next;
2490	}
2491	//slAct = slAct->next;
2492	//delete deleteMarker;
2493	}//if(!isParallel(fNormal, slNormal))
2494	if( (slAct->next==NULL) && (maybe==TRUE) )
2495	{
2496	doNotAdd=FALSE;
2497	break;
2498	}
2499	slAct = slAct->next;
2500	/* NOTE The following lines must not be here but rather called inside the if-block above.
2501	Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now,
2502	(not nice!) but since they are in seperate branches of the if-statement there should not
2503	be a way it gets called twice thus ommiting one facet:
2504	slAct = slAct->next;*/
2505	//delete deleteMarker;
2506	//deleteMarker=NULL;
2507	//if slAct was marked as to be deleted, delete it here!
2508	}//while(slAct!=NULL)
2509	if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) \|\| (doNotAdd==FALSE) )
2510	//if( (notParallelCtr==lengthOfSearchList ) \|\| doNotAdd==FALSE )
2511	{
2512	#ifdef gfan_DEBUG
2513	cout << "Adding facet (";
2514	fNormal->show(1,0);
2515	cout << ") to SLA " << endl;
2516	#endif
2517	//Add fAct to SLA
2518	facet *marker;
2519	marker = slEnd;
2520	facet *f2Act;
2521	f2Act = fAct->codim2Ptr;
2522
2523	slEnd->next = new facet();
2524	slEnd = slEnd->next;
2525	facet *slEndCodim2Root;
2526	facet *slEndCodim2Act;
2527	slEndCodim2Root = slEnd->codim2Ptr;
2528	slEndCodim2Act = slEndCodim2Root;
2529
2530	slEnd->setUCN(this->getUCN());
2531	slEnd->isFlippable = TRUE;
2532	slEnd->setFacetNormal(fNormal);
2533	slEnd->prev = marker;
2534	//Copy codim2-facets
2535	intvec *f2Normal;// = new intvec(this->numVars);
2536	while(f2Act!=NULL)
2537	{
2538	f2Normal=f2Act->getFacetNormal();
2539	if(slEndCodim2Root==NULL)
2540	{
2541	slEndCodim2Root = new facet(2);
2542	slEnd->codim2Ptr = slEndCodim2Root;
2543	slEndCodim2Root->setFacetNormal(f2Normal);
2544	slEndCodim2Act = slEndCodim2Root;
2545	}
2546	else
2547	{
2548	slEndCodim2Act->next = new facet(2);
2549	slEndCodim2Act = slEndCodim2Act->next;
2550	slEndCodim2Act->setFacetNormal(f2Normal);
2551	}
2552	f2Act = f2Act->next;
2553	}
2554	lengthOfSearchList++;
2555	}//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) \|\|
2556	fAct = fAct->next;
2557	}//if(fAct->isFlippable==TRUE)
2558	else
2559	{
2560	fAct = fAct->next;
2561	}
2562	}//while(fAct!=NULL)
2563	return slHead;
2564	// delete sl2Normal;
2565	// delete slNormal;
2566	// delete f2Normal;
2567	// delete fNormal;
2568	}//addC2N
2569
2570	/** \brief Compute the gcd of two ints
2571	*/
2572	int gcone::intgcd(int a, int b)
2573	{
2574	int r, p=a, q=b;
2575	if(p < 0)
2576	p = -p;
2577	if(q < 0)
2578	q = -q;
2579	while(q != 0)
2580	{
2581	r = p % q;
2582	p = q;
2583	q = r;
2584	}
2585	return p;
2586	}
2587
2588	/** \brief Construct a dd_MatrixPtr from a cone's list of facets
2589	*
2590	*/
2591	dd_MatrixPtr gcone::facets2Matrix(gcone const &gc)
2592	{
2593	facet *fAct;
2594	fAct = gc.facetPtr;
2595
2596	dd_MatrixPtr M;
2597	dd_rowrange ddrows;
2598	dd_colrange ddcols;
2599	ddcols=(this->numVars)+1;
2600	ddrows=this->numFacets;
2601	dd_NumberType numb = dd_Integer;
2602	M=dd_CreateMatrix(ddrows,ddcols);
2603
2604	int jj=0;
2605
2606	while(fAct!=NULL)
2607	{
2608	intvec *comp;
2609	comp = fAct->getFacetNormal();
2610	for(int ii=0;ii<this->numVars;ii++)
2611	{
2612	dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
2613	}
2614	jj++;
2615	fAct=fAct->next;
2616	}
2617
2618	return M;
2619	}
2620
2621	/** \brief Write information about a cone into a file on disk
2622	*
2623	* This methods writes the information needed for the "second" method into a file.
2624	* The file's is divided in sections containing information on
2625	* 1) the ring
2626	* 2) the cone's Groebner Basis
2627	* 3) the cone's facets
2628	* Each line contains exactly one date
2629	* Each section starts with its name in CAPITALS
2630	*/
2631	void gcone::writeConeToFile(gcone const &gc, bool usingIntPoints)
2632	{
2633	int UCN=gc.UCN;
2634	stringstream ss;
2635	ss << UCN;
2636	string UCNstr = ss.str();
2637
2638	string prefix="/tmp/cone";
2639	// string prefix="./"; //crude hack -> run tests in separate directories
2640	string suffix=".sg";
2641	string filename;
2642	filename.append(prefix);
2643	filename.append(UCNstr);
2644	filename.append(suffix);
2645
2646
2647	// int thisPID = getpid();
2648	// ss << thisPID;
2649	// string strPID = ss.str();
2650	// string prefix="./";
2651
2652	ofstream gcOutputFile(filename.c_str());
2653	facet *fAct;
2654	fAct = gc.facetPtr;
2655	facet *f2Act;
2656	f2Act=fAct->codim2Ptr;
2657
2658	char *ringString = rString(currRing);
2659
2660	if (!gcOutputFile)
2661	{
2662	cout << "Error opening file for writing in writeConeToFile" << endl;
2663	}
2664	else
2665	{
2666	gcOutputFile << "UCN" << endl;
2667	gcOutputFile << this->UCN << endl;
2668	gcOutputFile << "RING" << endl;
2669	gcOutputFile << ringString << endl;
2670	gcOutputFile << "GCBASISLENGTH" << endl;
2671	gcOutputFile << IDELEMS(gc.gcBasis) << endl;
2672	//Write this->gcBasis into file
2673	gcOutputFile << "GCBASIS" << endl;
2674	for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
2675	{
2676	gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
2677	}
2678
2679	gcOutputFile << "FACETS" << endl;
2680
2681	while(fAct!=NULL)
2682	{
2683	intvec *iv;
2684	intvec *iv2;
2685	iv=fAct->getFacetNormal();
2686	f2Act=fAct->codim2Ptr;
2687	for (int ii=0;ii<iv->length();ii++)
2688	{
2689	if (ii<iv->length()-1)
2690	{
2691	gcOutputFile << (*iv)[ii] << ",";
2692	}
2693	else
2694	{
2695	gcOutputFile << (*iv)[ii] << " ";
2696	}
2697	}
2698	while(f2Act!=NULL)
2699	{
2700	iv2=f2Act->getFacetNormal();
2701	for(int jj=0;jj<iv2->length();jj++)
2702	{
2703	if (jj<iv2->length()-1)
2704	{
2705	gcOutputFile << (*iv2)[jj] << ",";
2706	}
2707	else
2708	{
2709	gcOutputFile << (*iv2)[jj] << " ";
2710	}
2711	}
2712	f2Act = f2Act->next;
2713	}
2714	gcOutputFile << endl;
2715	fAct=fAct->next;
2716	}
2717	gcOutputFile.close();
2718	}
2719
2720	}//writeConeToFile(gcone const &gc)
2721
2722	/** \brief Reads a cone from a file identified by its number
2723	*/
2724	void gcone::readConeFromFile(int UCN, gcone *gc)
2725	{
2726	//int UCN=gc.UCN;
2727	stringstream ss;
2728	ss << UCN;
2729	string UCNstr = ss.str();
2730	string line;
2731	string strGcBasisLength;
2732	string strMonom, strCoeff, strCoeffNom, strCoeffDenom;
2733	int gcBasisLength=0;
2734	int intCoeff=1;
2735	int intCoeffNom=1; //better (number) but that's not compatible with stringstream;
2736	int intCoeffDenom=1;
2737	number nCoeff=nInit(1);
2738	number nCoeffNom=nInit(1);
2739	number nCoeffDenom=nInit(1);
2740	bool hasCoeffInQ = FALSE; //does the polynomial have rational coeff?
2741	bool hasNegCoeff = FALSE; //or a negative one?
2742	size_t found; //used for find_first_(not)_of
2743
2744	string prefix="/tmp/cone";
2745	string suffix=".sg";
2746	string filename;
2747	filename.append(prefix);
2748	filename.append(UCNstr);
2749	filename.append(suffix);
2750
2751	ifstream gcInputFile(filename.c_str(), ifstream::in);
2752
2753	ring saveRing=currRing;
2754	rChangeCurrRing(gc->baseRing);
2755	poly strPoly=pInit();
2756	poly resPoly=pInit(); //The poly to be read in
2757
2758	string::iterator EOL;
2759	int terms=1; //#Terms in the poly
2760
2761	while( !gcInputFile.eof() )
2762	{
2763	getline(gcInputFile,line);
2764	hasCoeffInQ = FALSE;
2765	hasNegCoeff = FALSE;
2766
2767	if(line=="GCBASISLENGTH")
2768	{
2769	getline(gcInputFile, line);
2770	strGcBasisLength = line;
2771	gcBasisLength=atoi(strGcBasisLength.c_str());
2772	}
2773	if(line=="GCBASIS")
2774	{
2775	for(int jj=0;jj<gcBasisLength;jj++)
2776	{
2777	getline(gcInputFile,line);
2778	//find out how many terms the poly consists of
2779	for(EOL=line.begin(); EOL!=line.end();++EOL)
2780	{
2781	string s;
2782	s=*EOL;
2783	if(s=="+" \|\| s=="-")
2784	terms++;
2785	}
2786	//magically convert strings into polynomials
2787	//polys.cc:p_Read
2788	//check until first occurance of + or -
2789	//data or c_str
2790	while(!line.empty())
2791	{
2792	hasNegCoeff = FALSE;
2793	found = line.find_first_of("+-"); //get the first monomial
2794	string tmp;
2795	tmp=line[found];
2796	// if(found!=0 && (tmp.compare("-")==0) )
2797	// hasNegCoeff = TRUE; //We have a coeff < 0
2798	if(found==0)
2799	{
2800	if(tmp.compare("-")==0)
2801	hasNegCoeff = TRUE;
2802	line.erase(0,1); //remove leading + or -
2803	found = line.find_first_of("+-"); //adjust found
2804	}
2805	strMonom = line.substr(0,found);
2806	line.erase(0,found);
2807	found = strMonom.find_first_of("/");
2808	if(found!=string::npos) //i.e. "/" exists in strMonom
2809	{
2810	hasCoeffInQ = TRUE;
2811	strCoeffNom=strMonom.substr(0,found);
2812	strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890"));
2813	strMonom.erase(0,found);
2814	strMonom.erase(0,strMonom.find_first_not_of("1234567890/"));
2815	// ss << strCoeffNom;
2816	// ss >> intCoeffNom;
2817	// nCoeffNom=(snumber*)strCoeffNom.c_str();
2818	nRead(strCoeffNom.c_str(), &nCoeffNom);
2819	nRead(strCoeffDenom.c_str(), &nCoeffDenom);
2820	// nCoeffNom=nInit(intCoeffNom);
2821	// ss << strCoeffDenom;
2822	// ss >> intCoeffDenom;
2823	// nCoeffDenom=nInit(intCoeffDenom);
2824	// nCoeffDenom=(snumber*)strCoeffNom.c_str();
2825	//NOTE NOT SURE WHETHER THIS WILL WORK!
2826	//nCoeffNom=nInit(intCoeffNom);
2827	// nCoeffDenom=(number)strCoeffDenom.c_str();
2828	// nCoeffDenom=(number)strCoeffDenom.c_str();
2829	}
2830	else
2831	{
2832	found = strMonom.find_first_not_of("1234567890");
2833	strCoeff = strMonom.substr(0,found);
2834	if(!strCoeff.empty())
2835	{
2836	nRead(strCoeff.c_str(),&nCoeff);
2837	}
2838	else
2839	{
2840	// intCoeff = 1;
2841	nCoeff = nInit(1);
2842	}
2843
2844	}
2845	const char* monom = strMonom.c_str();
2846
2847	p_Read(monom,strPoly,currRing);
2848	switch (hasCoeffInQ)
2849	{
2850	case TRUE:
2851	if(hasNegCoeff)
2852	nCoeffNom=nNeg(nCoeffNom);
2853	// intCoeffNom *= -1;
2854	// pSetCoeff(strPoly, nDiv((number)intCoeffNom, (number)intCoeffDenom));
2855	pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom));
2856	break;
2857	case FALSE:
2858	if(hasNegCoeff)
2859	nCoeff=nNeg(nCoeff);
2860	// intCoeff *= -1;
2861	if(!nIsOne(nCoeff))
2862	{
2863	// if(hasNegCoeff)
2864	// intCoeff *= -1;
2865	// pSetCoeff(strPoly,(number) intCoeff);
2866	pSetCoeff(strPoly, nCoeff );
2867	}
2868	break;
2869
2870	}
2871	//pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1???
2872	if(pIsConstantComp(resPoly))
2873	resPoly=pCopy(strPoly);
2874	else
2875	resPoly=pAdd(resPoly,strPoly);
2876
2877
2878	}//while(!line.empty())
2879
2880	gcBasis->m[jj]=pCopy(resPoly);
2881	resPoly=NULL; //reset
2882	}
2883	break;
2884	}
2885	}//while(!gcInputFile.eof())
2886
2887	gcInputFile.close();
2888	rChangeCurrRing(saveRing);
2889	}
2890
2891	// static void gcone::idPrint(ideal &I)
2892	// {
2893	// for(int ii=0;ii<IDELEMS(I);ii++)
2894	// {
2895	// cout << "_[" << ii << "]=" << I->m[ii] << endl;
2896	// }
2897	// }
2898
2899	void gcone::lprepareResult(lists l, gcone &gc)
2900	{
2901	gcone *gcAct;
2902	gcAct = &gc;
2903	matrix mFacetNormal=mpNew(counter,gc.numVars);
2904
2905	facet *fAct;
2906	fAct = gc.facetPtr;
2907	while( gcAct!=NULL )
2908	{
2909	l->m[0].data=(int*)gc.getUCN();
2910	l->m[1].data=(intvec*)(&gcAct->f2M(gcAct,gcAct->facetPtr));
2911	gcAct = gcAct->next;
2912	}
2913	}
2914	/** \brief Write facets of a cone into a matrix
2915	* Takes a pointer to a facet as 2nd arg in order to determine whether
2916	* it operates in codim 1 or 2
2917	*/
2918	intvec gcone::f2M(gcone gc, facet f)
2919	{
2920	facet *fAct;
2921	int codim=1;
2922	if(f==gc->facetPtr)
2923	fAct = gc->facetPtr;
2924	else
2925	{
2926	fAct = gc->facetPtr->codim2Ptr;
2927	codim=2;
2928	}
2929	intvec mRes=new intvec(counter,gc->numVars,0);//nrows==counter, ncols==numVars
2930	intvec *fNormal = new intvec(this->numVars);
2931	for(int ii=0;ii<gc->numFacets*(this->numVars);ii++)
2932	{
2933	fNormal=fAct->getFacetNormal();
2934	for(int jj=0;jj<this->numVars ;jj++ )
2935	{
2936	mRes[ii]=(*fNormal)[jj];
2937	ii++;
2938	}
2939	fAct = fAct->next;
2940	}
2941	return mRes;
2942	}
2943
2944	int gcone::counter=0;
2945	int gfanHeuristic;
2946	// ideal gfan(ideal inputIdeal, int h)
2947	lists gfan(ideal inputIdeal, int h)
2948	{
2949	lists lResList=(lists)omAllocBin(slists_bin);
2950	lResList->Init(5);
2951	lResList->m[0].rtyp=INT_CMD;
2952	lResList->m[0].data=(int*)255;
2953	lResList->m[1].rtyp=INTVEC_CMD;
2954	lResList->m[1].data=(intvec*)NULL;
2955	lResList->m[2].rtyp=IDEAL_CMD;
2956	lResList->m[2].data=(ideal*)NULL;
2957	lResList->m[3].rtyp=RING_CMD;
2958	lResList->m[3].data=(ring*)NULL;
2959	lResList->m[4].rtyp=LIST_CMD;
2960	lResList->m[4].data=(lists)NULL;
2961	int numvar = pVariables;
2962	gfanHeuristic = h;
2963
2964	enum searchMethod {
2965	reverseSearch,
2966	noRevS
2967	};
2968
2969	searchMethod method;
2970	method = noRevS;
2971	// method = reverseSearch;
2972
2973	#ifdef gfan_DEBUG
2974	cout << "Now in subroutine gfan" << endl;
2975	#endif
2976	ring inputRing=currRing; // The ring the user entered
2977	ring rootRing; // The ring associated to the target ordering
2978	ideal res;
2979	facet *fRoot;
2980
2981	if (method==reverseSearch)
2982	{
2983
2984	/* Construct a new ring which will serve as our root*/
2985	rootRing=rCopy0(currRing);
2986	rootRing->order[0]=ringorder_lp;
2987
2988	rComplete(rootRing);
2989	rChangeCurrRing(rootRing);
2990
2991	/* Fetch the inputIdeal into our rootRing */
2992	map theMap=(map)idMaxIdeal(1); //evil hack!
2993	theMap->preimage=NULL; //neccessary?
2994	ideal rootIdeal;
2995	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
2996	#ifndef NDEBUG
2997	#ifdef gfan_DEBUG
2998	cout << "Root ideal is " << endl;
2999	idShow(rootIdeal);
3000	cout << "The root ring is " << endl;
3001	rWrite(rootRing);
3002	cout << endl;
3003	#endif
3004	#endif
3005
3006	//gcone *gcRoot = new gcone(); //Instantiate the sink
3007	gcone *gcRoot = new gcone(rootRing,rootIdeal);
3008	gcone *gcAct;
3009	gcAct = gcRoot;
3010	gcAct->numVars=pVariables;
3011	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
3012	#ifndef NDEBUG
3013	idShow(gcAct->gcBasis);
3014	#endif
3015	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
3016	//gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
3017	/*Now it is time to compute the search facets, respectively start the reverse search.
3018	But since we are in the root all facets should be search facets. IS THIS TRUE?
3019	NOTE: Check for flippability is not very sophisticated
3020	*/
3021	//gcAct->reverseSearch(gcAct);
3022	rChangeCurrRing(rootRing);
3023	res=gcRoot->gcBasis;
3024	}//if method==reverSearch
3025
3026	if(method==noRevS)
3027	{
3028	gcone *gcRoot = new gcone(currRing,inputIdeal);
3029	gcone *gcAct;
3030	gcAct = gcRoot;
3031	gcAct->numVars=pVariables;
3032	gcAct->getGB(inputIdeal);
3033	cout << "GB of input ideal is:" << endl;
3034	#ifndef NDEBUG
3035	idShow(gcAct->gcBasis);
3036	#endif
3037	if(gcAct->isMonomial(gcAct->gcBasis))
3038	{
3039	WerrorS("Monomial input - terminating");
3040	res=gcAct->gcBasis;
3041	//break;
3042	}
3043	cout << endl;
3044	gcAct->getConeNormals(gcAct->gcBasis);
3045	gcAct->noRevS(*gcAct);
3046	//res=gcAct->gcBasis;
3047	//Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
3048	res = inputIdeal;
3049	// intvec mRes=gcRoot->f2M(gcRoot,gcRoot->facetPtr);
3050	gcRoot->lprepareResult(lResList,*gcRoot);
3051	}
3052	dd_free_global_constants();
3053
3054	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
3055	The return type will then be of type LIST_CMD
3056	Assume gfan has finished, thus we have enumerated all the cones
3057	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
3058	Groebner Basis and merge this somehow with LIST_CMD
3059	=> Count the cones!
3060	*/
3061	//rChangeCurrRing(rootRing);
3062	//res=gcAct->gcBasis;
3063	//res=gcRoot->gcBasis;
3064	// return res;
3065	return lResList;
3066	//return GBlist;
3067	}
3068	/*
3069	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
3070	*/
3071	#endif

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

Download in other formats: