Context Navigation

source: git/kernel/gfan.cc @ a7efc97

spielwiese

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

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

Download in other formats: