Context Navigation

source: git/kernel/gfan.cc @ e7bc2f

spielwiese

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

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

Download in other formats: