Context Navigation

source: git/kernel/gfan.cc @ 7c0ec6

spielwiese

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

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

Download in other formats: