Context Navigation

source: git/kernel/gfan.cc @ 11a7dc

spielwiese

Last change on this file since 11a7dc was 11a7dc, checked in by Martin Monerjan, 14 years ago
sed s/cout/printf/g Removed facets2Matrix and interiorPoint2 git-svn-id: file:///usr/local/Singular/svn/trunk@12649 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 122.2 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	#include "options.h"
13	#include "kstd1.h"
14	#include "kutil.h" //ksCreateSpoly
15	// #include "int64vec.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 "stairc.h" //For Hilbert series
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 <stdlib.h>
34	//#include <gmpxx.h>
35	// #include <vector>
36	#include <assert.h>
37
38
39	/DO NOT REMOVE THIS/
40	#ifndef GMPRATIONAL
41	#define GMPRATIONAL
42	#endif
43
44	//Hacks for different working places
45	#define p800
46
47	#ifdef UNI
48	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
49	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
50	#endif
51
52	#ifdef HOME
53	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
54	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
55	#endif
56
57	#ifdef ITWM
58	#include "/u/slg/monerjan/cddlib/include/setoper.h"
59	#include "/u/slg/monerjan/cddlib/include/cdd.h"
60	#include "/u/slg/monerjan/cddlib/include/cddmp.h"
61	#endif
62
63	#ifdef p800
64	#include "../../cddlib/include/setoper.h"
65	#include "../../cddlib/include/cdd.h"
66	#include "../../cddlib/include/cddmp.h"
67	#endif
68
69	#ifndef gfan_DEBUG
70	// #define gfan_DEBUG
71	#ifndef gfan_DEBUGLEVEL
72	#define gfan_DEBUGLEVEL 1
73	#endif
74	#endif
75
76	//NOTE Defining this will slow things down!
77	//Only good for very coarse profiling
78	#define gfanp
79	#ifdef gfanp
80	#include <sys/time.h>
81	#include <iostream>
82	#endif
83
84	#ifndef SHALLOW
85	#define SHALLOW
86	#endif
87
88	#include <gfan.h>
89	using namespace std;
90
91	#define ivIsStrictlyPositive iv64isStrictlyPositive
92
93	/**
94	*\brief Class facet
95	* Implements the facet structure as a linked list
96	*
97	*/
98
99	/** \brief The default constructor for facets
100	*/
101	facet::facet()
102	{
103	// Pointer to next facet. */
104	/* Defaults to NULL. This way there is no need to check explicitly */
105	this->fNormal=NULL;
106	this->interiorPoint=NULL;
107	this->UCN=0;
108	this->codim2Ptr=NULL;
109	this->codim=1; //default for (codim-1)-facets
110	this->numCodim2Facets=0;
111	this->numRays=0;
112	this->flipGB=NULL;
113	this->next=NULL;
114	this->prev=NULL;
115	this->flipRing=NULL; //the ring on the other side
116	this->isFlippable=FALSE;
117	}
118
119	/** \brief Constructor for facets of codim == 2
120	* Note that as of now the code of the constructors is only for facets and codim2-faces. One
121	* could easily change that by renaming numCodim2Facets to numCodimNminusOneFacets or similar
122	*/
123	facet::facet(const int &n)
124	{
125	this->fNormal=NULL;
126	this->interiorPoint=NULL;
127	this->UCN=0;
128	this->codim2Ptr=NULL;
129	if(n==2)
130	{
131	this->codim=n;
132	}//NOTE Handle exception here!
133	this->numCodim2Facets=0;
134	this->numRays=0;
135	this->flipGB=NULL;
136	this->next=NULL;
137	this->prev=NULL;
138	this->flipRing=NULL;
139	this->isFlippable=FALSE;
140	}
141
142	/** \brief The copy constructor
143	* By default only copies the fNormal, f2Normals and UCN
144	*/
145	facet::facet(const facet& f)
146	{
147	this->fNormal=ivCopy(f.fNormal);
148	this->UCN=f.UCN;
149	this->isFlippable=f.isFlippable;
150	//Needed for flip2
151	//NOTE ONLY REFERENCE
152	this->interiorPoint=ivCopy(f.interiorPoint);//only referencing is prolly not the best thing to do in a copy constructor
153	facet* f2Copy;
154	f2Copy=f.codim2Ptr;
155	facet* f2Act;
156	f2Act=this->codim2Ptr;
157	while(f2Copy!=NULL)
158	{
159	if(f2Act==NULL \|\| f2Act==(facet*)0xfefefefefefefefe)
160	{
161	f2Act=new facet(2);
162	this->codim2Ptr=f2Act;
163	}
164	else
165	{
166	facet* marker;
167	marker = f2Act;
168	f2Act->next = new facet(2);
169	f2Act = f2Act->next;
170	f2Act->prev = marker;
171	}
172	intvec *f2Normal;
173	f2Normal = f2Copy->getFacetNormal();
174	// f2Act->setFacetNormal(f2Copy->getFacetNormal());
175	f2Act->setFacetNormal(f2Normal);
176	delete f2Normal;
177	f2Act->setUCN(f2Copy->getUCN());
178	f2Copy = f2Copy->next;
179	}
180	}
181
182	/** \brief Shallow copy constructor for facets
183	* We only need the interior point for equality testing
184	*/
185	facet* facet::shallowCopy(const facet& f)
186	{
187	facet *res = new facet();
188	res->fNormal=(intvec * const)f.fNormal;
189	res->UCN=f.UCN;
190	res->isFlippable=f.isFlippable;
191	res->interiorPoint=(intvec * const)f.interiorPoint;
192	res->codim2Ptr=(facet * const)f.codim2Ptr;
193	res->prev=NULL;
194	res->next=NULL;
195	res->flipGB=NULL;
196	res->flipRing=NULL;
197	return res;
198	}
199
200	void facet::shallowDelete()
201	{
202	#ifdef gfan_DEBUG
203	printf("shallowdel@UCN %i\n", this->getUCN());
204	#endif
205	this->fNormal=NULL;
206	// this->UCN=0;
207	this->interiorPoint=NULL;
208	this->codim2Ptr=NULL;
209	this->prev=NULL;
210	this->next=NULL;
211	this->flipGB=NULL;
212	this->flipRing=NULL;
213	assert(this->fNormal==NULL);
214	// delete(this);
215	}
216
217	/** The default destructor */
218	facet::~facet()
219	{
220	#ifdef gfan_DEBUG
221	printf("~facet@UCN %i\n",this->getUCN());
222	#endif
223	if(this->fNormal!=NULL)
224	delete this->fNormal;
225	if(this->interiorPoint!=NULL)
226	delete this->interiorPoint;
227	/* Cleanup the codim2-structure */
228	if(this->codim==2)
229	// if(this->codim2Ptr!=NULL)
230	{
231	facet *codim2Ptr;
232	codim2Ptr = this->codim2Ptr;
233	while(codim2Ptr!=NULL)
234	{
235	if(codim2Ptr->fNormal!=NULL)
236	{
237	delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
238	codim2Ptr = codim2Ptr->next;
239	}
240	}
241	// delete this->codim2Ptr;
242	}
243	//The rays are stored in the cone!
244	// if(this->codim2Ptr!=NULL)
245	// delete this->codim2Ptr;
246	if(this->flipGB!=NULL)
247	idDelete((ideal *)&this->flipGB);
248	if(this->flipRing!=NULL && this->flipRing->idroot!=(idhdl)0xfbfbfbfbfbfbfbfb)
249	// rDelete(this->flipRing); //See vol II/134
250	// this->flipRing=NULL;
251	this->prev=NULL;
252	this->next=NULL;
253	}
254
255	inline const intvec *facet::getRef2FacetNormal() const
256	{
257	return(this->fNormal);
258	}
259
260	/** Equality check for facets based on unique interior points*/
261	static bool areEqual2(facet* f, facet *g)
262	{
263	#ifdef gfanp
264	gcone::numberOfFacetChecks++;
265	timeval start, end;
266	gettimeofday(&start, 0);
267	#endif
268	bool res = TRUE;
269	const intvec *fIntP = f->getRef2InteriorPoint();
270	const intvec *gIntP = g->getRef2InteriorPoint();
271	for(int ii=0;ii<pVariables;ii++)
272	{
273	if( (fIntP)[ii] != (gIntP)[ii] )
274	{
275	res=FALSE;
276	break;
277	}
278	}
279	#ifdef gfanp
280	gettimeofday(&end, 0);
281	gcone::t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
282	#endif
283	return res;
284	}
285
286	/** \brief Comparison of facets
287	* called from enqueueNewFacets
288	* The facet normals are primitve vectors since we call gcone::normalize() on all cones.
289	* Hence it should suffice to check whether facet normal f equals minus facet normal s.
290	* If so we check the extremal rays
291	*
292	* BEWARE: It would be better to use const intvec* but that will lead to call something like
293	* int foo=((intvec)f2Normal)->compare((intvec)s2Normal) resulting in much higher memory usage
294	*/
295	static bool areEqual(facet f, facet s)
296	{
297	#ifdef gfanp
298	gcone::numberOfFacetChecks++;
299	timeval start, end;
300	gettimeofday(&start, 0);
301	#endif
302	bool res = TRUE;
303	int notParallelCtr=0;
304	int ctr=0;
305	const intvec* fNormal; //No new since ivCopy and therefore getFacetNormal return a new
306	const intvec* sNormal;
307	fNormal = f->getRef2FacetNormal();//->getFacetNormal();
308	sNormal = s->getRef2FacetNormal();//->getFacetNormal();
309	//Do not need parallelity. Too time consuming
310	// if(!isParallel(fNormal,sNormal))
311	// if(fNormal->compare(ivNeg(sNormal))!=0)//This results in a Mandelbug
312	// notParallelCtr++;
313	// else//parallelity, so we check the codim2-facets
314	intvec fNRef=const_cast<intvec>(fNormal);
315	intvec sNRef=const_cast<intvec>(sNormal);
316	// if(isParallel(fNormal,sNormal))
317	if(isParallel(fNRef,sNRef))
318	// if(fNormal->compare((sNormal))!=0)//Behold! Teh definitive Mandelbug
319	{
320	facet* f2Act;
321	facet* s2Act;
322	f2Act = f->codim2Ptr;
323	ctr=0;
324	while(f2Act!=NULL)
325	{
326	const intvec* f2Normal;
327	f2Normal = f2Act->getRef2FacetNormal();//->getFacetNormal();
328	// intvec f2Ref=const_cast<intvec>(f2Normal);
329	s2Act = s->codim2Ptr;
330	while(s2Act!=NULL)
331	{
332	const intvec* s2Normal;
333	s2Normal = s2Act->getRef2FacetNormal();//->getFacetNormal();
334	// bool foo=areEqual(f2Normal,s2Normal);
335	// intvec s2Ref=const_cast<intvec>(s2Normal);
336	int foo=f2Normal->compare(s2Normal);
337	if(foo==0)
338	ctr++;
339	s2Act = s2Act->next;
340	// delete s2Normal;
341	}
342	// delete f2Normal;
343	f2Act = f2Act->next;
344	}
345	}
346	// delete fNormal;
347	// delete sNormal;
348	if(ctr==f->numCodim2Facets)
349	res=TRUE;
350	else
351	{
352	#ifdef gfanp
353	gcone::parallelButNotEqual++;
354	#endif
355	res=FALSE;
356	}
357	#ifdef gfanp
358	gettimeofday(&end, 0);
359	gcone::t_areEqual += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
360	#endif
361	return res;
362	// intvec *foo=ivNeg(sNormal);
363	// if(fNormal->compare(foo)!=0) //facet normals
364	// {
365	// delete foo;
366	// res=FALSE;
367	// }
368	// else
369	// {
370	// facet* f2Act;
371	// facet* s2Act;
372	// f2Act = f->codim2Ptr;
373	// ctr=0;
374	// while(f2Act!=NULL)
375	// {
376	// intvec* f2Normal;
377	// f2Normal = f2Act->getFacetNormal();
378	// s2Act = s->codim2Ptr;
379	// while(s2Act!=NULL)
380	// {
381	// intvec* s2Normal;
382	// s2Normal = s2Act->getFacetNormal();
383	// // bool foo=areEqual(f2Normal,s2Normal);
384	// int foo=f2Normal->compare(s2Normal);
385	// if(foo==0)
386	// ctr++;
387	// s2Act = s2Act->next;
388	// delete s2Normal;
389	// }
390	// delete f2Normal;
391	// f2Act = f2Act->next;
392	// }
393	// }
394	// delete fNormal;
395	// delete sNormal;
396	// if(ctr==f->numCodim2Facets)
397	// res=TRUE;
398	// return res;
399	}
400
401	/** Stores the facet normal \param intvec*/
402	inline void facet::setFacetNormal(intvec *iv)
403	{
404	if(this->fNormal!=NULL)
405	delete this->fNormal;
406	this->fNormal = ivCopy(iv);
407	}
408
409	/** Hopefully returns the facet normal
410	* Mind: ivCopy returns a new intvec, so use this in the following way:
411	* intvec *iv;
412	* iv = this->getFacetNormal();
413	* [...]
414	* delete(iv);
415	*/
416	inline intvec *facet::getFacetNormal() const
417	{
418	return ivCopy(this->fNormal);
419	// return this->fNormal;
420	}
421
422	/** Method to print the facet normal*/
423	inline void facet::printNormal() const
424	{
425	fNormal->show();
426	}
427
428	/** Store the flipped GB*/
429	inline void facet::setFlipGB(ideal I)
430	{
431	this->flipGB=idCopy(I);
432	}
433
434	/** Returns a pointer to the flipped GB
435	Seems not be used
436	Anyhow idCopy would make sense here.
437	*/
438	inline ideal facet::getFlipGB()
439	{
440	return this->flipGB;
441	}
442
443	/** Print the flipped GB*/
444	inline void facet::printFlipGB()
445	{
446	#ifndef NDEBUG
447	idShow(this->flipGB);
448	#endif
449	}
450
451	/** Set the UCN */
452	inline void facet::setUCN(int n)
453	{
454	this->UCN=n;
455	}
456
457	/** \brief Get the UCN
458	* Returns the UCN iff this != NULL, else -1
459	*/
460	inline int facet::getUCN()
461	{
462	#ifndef NDEBUG
463	#if SIZEOF_LONG==8
464	if((this!=NULL && this!=(facet * const)0xfbfbfbfbfbfbfbfb))
465	#elif SIZEOF_LONG==4
466	if((this!=NULL && this!=(facet * const)0xfbfbfbfb))
467	#endif
468	#endif
469	#ifdef NDEBUG
470	if(this!=NULL)
471	#endif
472	return this->UCN;
473	else
474	return -1;
475	}
476
477	/** Store an interior point of the facet */
478	inline void facet::setInteriorPoint(intvec *iv)
479	{
480	if(this->interiorPoint!=NULL)
481	delete this->interiorPoint;
482	this->interiorPoint = ivCopy(iv);
483	}
484
485	/** Returns a pointer to this->interiorPoint
486	* MIND: ivCopy returns a new intvec
487	* @see facet::getFacetNormal
488	*/
489	inline intvec *facet::getInteriorPoint()
490	{
491	return ivCopy(this->interiorPoint);
492	}
493
494	inline const intvec *facet::getRef2InteriorPoint()
495	{
496	return (this->interiorPoint);
497	}
498
499	/** \brief Debugging function
500	* prints the facet normal an all (codim-2)-facets that belong to it
501	*/
502	volatile void facet::fDebugPrint()
503	{
504	facet *codim2Act;
505	codim2Act = this->codim2Ptr;
506	intvec *fNormal;
507	fNormal = this->getFacetNormal();
508	printf("=======================\n");
509	printf("Facet normal = (");fNormal->show(1,1);printf(")\n");
510	printf("-----------------------\n");
511	printf("Codim2 facets:\n");
512	while(codim2Act!=NULL)
513	{
514	intvec *f2Normal;
515	f2Normal = codim2Act->getFacetNormal();
516	printf("(");f2Normal->show(1,0);printf(")\n");
517	codim2Act = codim2Act->next;
518	delete f2Normal;
519	}
520	printf("=======================\n");
521	delete fNormal;
522	}
523
524	//friend class gcone; //Bad style
525
526
527	/**
528	*\brief Implements the cone structure
529	*
530	* A cone is represented by a linked list of facet normals
531	* @see facet
532	*/
533
534
535	/** \brief Default constructor.
536	*
537	* Initialises this->next=NULL and this->facetPtr=NULL
538	*/
539	gcone::gcone()
540	{
541	this->next=NULL;
542	this->prev=NULL;
543	this->facetPtr=NULL; //maybe this->facetPtr = new facet();
544	this->baseRing=currRing;
545	this->counter++;
546	this->UCN=this->counter;
547	this->numFacets=0;
548	this->ivIntPt=NULL;
549	this->gcRays=NULL;
550	}
551
552	/** \brief Constructor with ring and ideal
553	*
554	* This constructor takes the root ring and the root ideal as parameters and stores
555	* them in the private members gcone::rootRing and gcone::inputIdeal
556	* This constructor is only called once in the computation of the GrÃ¶bner fan,
557	* namely for the very first cone. Therefore pred is set to 1.
558	* Might set it to this->UCN though...
559	* Since knowledge of the root ring is only needed when using reverse search,
560	* this constructor is not needed when using the "second" method
561	*/
562	gcone::gcone(ring r, ideal I)
563	{
564	this->next=NULL;
565	this->prev=NULL;
566	this->facetPtr=NULL;
567	this->inputIdeal=I;
568	this->baseRing=currRing;
569	this->counter++;
570	this->UCN=this->counter;
571	this->pred=1;
572	this->numFacets=0;
573	this->ivIntPt=NULL;
574	this->gcRays=NULL;
575	}
576
577	/** \brief Copy constructor
578	*
579	* Copies a cone, sets this->gcBasis to the flipped GB
580	* Call this only after a successful call to gcone::flip which sets facet::flipGB
581	*/
582	gcone::gcone(const gcone& gc, const facet &f)
583	{
584	this->next=NULL;
585	// this->prev=(gcone *)&gc; //comment in to get a tree
586	this->prev=NULL;
587	this->numVars=gc.numVars;
588	this->counter++;
589	this->UCN=this->counter;
590	this->pred=gc.UCN;
591	this->facetPtr=NULL;
592	this->gcBasis=idrCopyR(f.flipGB, f.flipRing);
593	// this->inputIdeal=idCopy(this->gcBasis);
594	this->baseRing=rCopy(f.flipRing);
595	this->numFacets=0;
596	this->ivIntPt=NULL;
597	this->gcRays=NULL;
598	}
599
600	/** \brief Default destructor
601	*/
602	gcone::~gcone()
603	{
604	#ifndef NDEBUG
605	#if SIZEOF_LONG==8
606	if( ( this->gcBasis!=(ideal)(0xfbfbfbfbfbfbfbfb) ) && (this->gcBasis!=NULL) )
607	idDelete((ideal*)&this->gcBasis);
608	#elif SIZEOF_LONG!=8
609	if(this->gcBasis!=(ideal)0xfbfbfbfb)
610	idDelete((ideal *)&this->gcBasis);
611	#endif
612	#else
613	if(this->gcBasis!=NULL)
614	idDelete((ideal *)&this->gcBasis);
615	#endif
616	// idDelete((ideal *)&this->gcBasis);
617	// if(this->inputIdeal!=NULL)
618	// idDelete((ideal *)&this->inputIdeal);
619	// if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
620	// rDelete(this->rootRing);
621	if(this->UCN!=1 && this->baseRing!=NULL)
622	rDelete(this->baseRing);
623	facet *fAct;
624	facet *fDel;
625	/Delete the facet structure/
626	fAct=this->facetPtr;
627	fDel=fAct;
628	while(fAct!=NULL)
629	{
630	fDel=fAct;
631	fAct=fAct->next;
632	delete fDel;
633	}
634	this->counter--;
635	//should be deleted in noRevS
636	// dd_FreeMatrix(this->ddFacets);
637	//dd_FreeMatrix(this->ddFacets);
638	for(int ii=0;ii<this->numRays;ii++)
639	delete(gcRays[ii]);
640	omFree(gcRays);
641	}
642
643	/** Returns the number of cones existing at the time*/
644	inline int gcone::getCounter()
645	{
646	return this->counter;
647	}
648
649	/** \brief Set the interior point of a cone */
650	inline void gcone::setIntPoint(intvec *iv)
651	{
652	if(this->ivIntPt!=NULL)
653	delete this->ivIntPt;
654	this->ivIntPt=ivCopy(iv);
655	}
656
657	/** \brief Returns either a physical copy the interior point of a cone or just a reference to it.*/
658	inline intvec *gcone::getIntPoint(bool shallow)
659	{
660	if(shallow==TRUE)
661	return this->ivIntPt;
662	else
663	return ivCopy(this->ivIntPt);
664	}
665
666	/** \brief Print the interior point */
667	inline void gcone::showIntPoint()
668	{
669	ivIntPt->show();
670	}
671
672	/** \brief Print facets
673	* This is mainly for debugging purposes. Usually called from within gdb
674	*/
675	volatile void gcone::showFacets(const short codim)
676	{
677	facet *f=this->facetPtr;
678	facet *f2=NULL;
679	if(codim==2)
680	f2=this->facetPtr->codim2Ptr;
681	while(f!=NULL)
682	{
683	intvec *iv;
684	iv = f->getFacetNormal();
685	printf("(");iv->show(1,0);
686	if(f->isFlippable==FALSE)
687	printf(")* ");
688	else
689	printf(") ");
690	delete iv;
691	if(codim==2)
692	{
693	f2=f->codim2Ptr;
694	while(f2!=NULL)
695	{
696	printf("[");f2->getFacetNormal()->show(1,0);printf("]");
697	f2 = f2->next;
698	}
699	printf("\n");
700	}
701	f=f->next;
702	}
703	printf("\n");
704	}
705
706	/** For debugging purposes only */
707	static volatile void showSLA(facet &f)
708	{
709	facet *fAct;
710	fAct = &f;
711	if(fAct!=NULL)
712	{
713	facet *codim2Act;
714	codim2Act = fAct->codim2Ptr;
715
716	printf("\n");
717	while(fAct!=NULL)
718	{
719	intvec *fNormal;
720	fNormal=fAct->getFacetNormal();
721	printf("(");fNormal->show(1,0);
722	if(fAct->isFlippable==TRUE)
723	printf(") ");
724	else
725	printf(")* ");
726	delete fNormal;
727	codim2Act = fAct->codim2Ptr;
728	printf(" Codim2: ");
729	while(codim2Act!=NULL)
730	{
731	intvec *f2Normal;
732	f2Normal = codim2Act->getFacetNormal();
733	printf("(");f2Normal->show(1,0);printf(") ");
734	delete f2Normal;
735	codim2Act = codim2Act->next;
736	}
737	printf("UCN = %i",fAct->getUCN());
738	fAct = fAct->next;
739	}
740	}
741	}
742
743	static void idDebugPrint(const ideal &I)
744	{
745	int numElts=IDELEMS(I);
746	printf("Ideal with %i generators\n", numElts);
747	printf("Leading terms: ");
748	for (int ii=0;ii<numElts;ii++)
749	{
750	pWrite0(pHead(I->m[ii]));
751	printf(",");
752	}
753	printf("\n");
754	}
755
756	static void invPrint(const ideal &I)
757	{
758	// int numElts=IDELEMS(I);
759	// cout << "inv = ";
760	// for(int ii=0;ii<numElts;ii++);
761	// {
762	// pWrite0(pHead(I->m[ii]));
763	// cout << ",";
764	// }
765	// cout << endl;
766	}
767
768	static bool isMonomial(const ideal &I)
769	{
770	bool res = TRUE;
771	for(int ii=0;ii<IDELEMS(I);ii++)
772	{
773	if(pLength((poly)I->m[ii])>1)
774	{
775	res = FALSE;
776	break;
777	}
778	}
779	return res;
780	}
781
782	/** \brief Set gcone::numFacets */
783	inline void gcone::setNumFacets()
784	{
785	}
786
787	/** \brief Get gcone::numFacets */
788	inline int gcone::getNumFacets()
789	{
790	return this->numFacets;
791	}
792
793	inline int gcone::getUCN()
794	{
795	if( this!=NULL)// && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
796	return this->UCN;
797	else
798	return -1;
799	}
800
801	inline int gcone::getPredUCN()
802	{
803	return this->pred;
804	}
805	/** Returns a copy of the this->baseRing */
806	inline ring gcone::getBaseRing()
807	{
808	return rCopy(this->baseRing);
809	}
810
811	inline ring gcone::getRef2BaseRing()
812	{
813	return this->baseRing;
814	}
815
816	/** \brief Compute the normals of the cone
817	*
818	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
819	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
820	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
821	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
822	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
823	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
824	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
825	*
826	* Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
827	* an interior point of the cone.
828	*/
829	void gcone::getConeNormals(const ideal &I, bool compIntPoint)
830	{
831	#ifdef gfanp
832	timeval start, end;
833	gettimeofday(&start, 0);
834	#endif
835	poly aktpoly;
836	int rows; // will contain the dimensions of the ineq matrix - deprecated by
837	dd_rowrange ddrows;
838	dd_colrange ddcols;
839	dd_rowset ddredrows; // # of redundant rows in ddineq
840	dd_rowset ddlinset; // the opposite
841	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
842	dd_NumberType ddnumb=dd_Integer; //Number type
843	dd_ErrorType dderr=dd_NoError;
844	//Compute the # inequalities i.e. rows of the matrix
845	rows=0; //Initialization
846	for (int ii=0;ii<IDELEMS(I);ii++)
847	{
848	// aktpoly=(poly)I->m[ii];
849	// rows=rows+pLength(aktpoly)-1;
850	rows=rows+pLength((poly)I->m[ii])-1;
851	}
852
853	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
854	ddrows=rows;
855	ddcols=this->numVars;
856	dd_MatrixPtr ddineq; //Matrix to store the inequalities
857	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
858
859	// We loop through each g\in GB and compute the resulting inequalities
860	for (int i=0; i<IDELEMS(I); i++)
861	{
862	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
863	//simpler version of storing expvect diffs
864	int leadexpv=(int)omAlloc(((this->numVars)+1)*sizeof(int));
865	// int tailexpv=(int)omAlloc(((this->numVars)+1)*sizeof(int));
866	pGetExpV(aktpoly,leadexpv);
867	while(pNext(aktpoly)!=NULL)
868	{
869	aktpoly=pNext(aktpoly);
870	int tailexpv=(int)omAlloc(((this->numVars)+1)*sizeof(int));
871	pGetExpV(aktpoly,tailexpv);
872	for(int kk=1;kk<=this->numVars;kk++)
873	{
874	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk],leadexpv[kk]-tailexpv[kk]);
875	}
876	aktmatrixrow += 1;
877	omFree(tailexpv);
878	}
879	omFree(leadexpv);
880	} //for
881	#if true
882	/*Let's make the preprocessing here. This could already be done in the above for-loop,
883	* but for a start it is more convenient here.
884	* We check the necessary condition of FJT p.18
885	* Quote: [...] every non-zero spoly should have at least one of its terms in inv(G)
886	*/
887	// ideal initialForm=idInit(IDELEMS(I),1);
888	// intvec *gamma=new intvec(this->numVars);
889	int falseGammaCounter=0;
890	int *redRowsArray=NULL;
891	int num_alloc=0;
892	int num_elts=0;
893	for(int ii=0;ii<ddineq->rowsize;ii++)
894	{
895	ideal initialForm=idInit(IDELEMS(I),1);
896	//read row ii into gamma
897	double tmp;
898	intvec *gamma=new intvec(this->numVars);
899	for(int jj=1;jj<=this->numVars;jj++)
900	{
901	tmp=mpq_get_d(ddineq->matrix[ii][jj]);
902	(*gamma)[jj-1]=(int)tmp;
903	}
904	computeInv((ideal&)I,initialForm,*gamma);
905	//Create leading ideal
906	ideal L=idInit(IDELEMS(initialForm),1);
907	for(int jj=0;jj<IDELEMS(initialForm);jj++)
908	{
909	L->m[jj]=pCopy(pHead(initialForm->m[jj]));
910	}
911
912	LObject *P = new sLObject();
913	memset(P,0,sizeof(LObject));
914
915	for(int jj=0;jj<=IDELEMS(initialForm)-2;jj++)
916	{
917	bool isMaybeFacet=FALSE;
918	P->p1=initialForm->m[jj]; //build spolys of initialForm in_v
919
920	for(int kk=jj+1;kk<=IDELEMS(initialForm)-1;kk++)
921	{
922	P->p2=initialForm->m[kk];
923	ksCreateSpoly(P);
924	if(P->p!=NULL) //spoly non zero=?
925	{
926	poly p=pInit();
927	poly q=pInit();
928	p=pCopy(P->p);
929	q=pHead(p); //Monomial q
930	isMaybeFacet=FALSE;
931	//TODO: Suffices to check LTs here
932	while(p!=NULL)
933	{
934	q=pHead(p);
935	// unsigned long sevSpoly=pGetShortExpVector(q);
936	// unsigned long not_sevL;
937	for(int ll=0;ll<IDELEMS(L);ll++)
938	{
939	// not_sevL=~pGetShortExpVector(L->m[ll]);//
940	//if(!(sevSpoly & not_sevL) && pLmDivisibleBy(L->m[ll],q) )//i.e. spoly is in L
941	if(pLmEqual(L->m[ll],q) \|\| pDivisibleBy(L->m[ll],q))
942	{
943	isMaybeFacet=TRUE;
944	break;//for
945	}
946	}
947	if(isMaybeFacet==TRUE)
948	{
949	break;//while(p!=NULL)
950	}
951	p=pNext(p);
952	}//while
953	pDelete(&p);
954	pDelete(&q);
955	if(isMaybeFacet==FALSE)
956	{
957	dd_set_si(ddineq->matrix[ii][0],1);
958	if(num_alloc==0)
959	num_alloc += 1;
960	else
961	num_alloc += 1;
962	void tmp = realloc(redRowsArray,(num_allocsizeof(int)));
963	if(!tmp)
964	{
965	WerrorS("Woah dude! Couldn't realloc memory\n");
966	exit(-1);
967	}
968	redRowsArray = (int*)tmp;
969	redRowsArray[num_elts]=ii;
970	num_elts++;
971	//break;//for(int kk, since we have found one that is not in L
972	goto _start; //mea culpa, mea culpa, mea maxima culpa
973	}
974	}//if(P->p!=NULL)
975	}//for k
976	}//for jj
977	_start:;
978	idDelete(&L);
979	delete P;
980	idDelete(&initialForm);
981	//idDelete(L);
982	delete gamma;
983	}//for(ii<ddineq-rowsize
984	// delete gamma;
985	int offset=0;//needed for correction of redRowsArray[ii]
986	for( int ii=0;ii<num_elts;ii++ )
987	{
988	dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
989	offset++;
990	}
991	free(redRowsArray);
992	/*And now for the strictly positive rows
993	* Doesn't gain significant speedup
994	*/
995	/int posRowsArray=NULL;
996	num_alloc=0;
997	num_elts=0;
998	for(int ii=0;ii<ddineq->rowsize;ii++)
999	{
1000	intvec *ivPos = new intvec(this->numVars);
1001	for(int jj=0;jj<this->numVars;jj++)
1002	(*ivPos)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
1003	bool isStrictlyPos=FALSE;
1004	int posCtr=0;
1005	for(int jj=0;jj<this->numVars;jj++)
1006	{
1007	intvec *ivCanonical = new intvec(this->numVars);
1008	jj==0 ? (ivCanonical)[ivPos->length()-1]=1 : (ivCanonical)[jj-1]=1;
1009	if(dotProduct(ivCanonical,ivPos)!=0)
1010	{
1011	if ((*ivPos)[jj]>=0)
1012	{
1013	posCtr++;
1014	}
1015	}
1016	delete ivCanonical;
1017	}
1018	if(posCtr==ivPos->length())
1019	isStrictlyPos=TRUE;
1020	if(isStrictlyPos==TRUE)
1021	{
1022	if(num_alloc==0)
1023	num_alloc += 1;
1024	else
1025	num_alloc += 1;
1026	void tmp = realloc(posRowsArray,(num_allocsizeof(int)));
1027	if(!tmp)
1028	{
1029	WerrorS("Woah dude! Couldn't realloc memory\n");
1030	exit(-1);
1031	}
1032	posRowsArray = (int*)tmp;
1033	posRowsArray[num_elts]=ii;
1034	num_elts++;
1035	}
1036	delete ivPos;
1037	}
1038	offset=0;
1039	for(int ii=0;ii<num_elts;ii++)
1040	{
1041	dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
1042	offset++;
1043	}
1044	free(posRowsArray);*/
1045	#endif
1046
1047	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
1048	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
1049	ddcols = ddineq->colsize;
1050
1051	this->ddFacets = dd_CopyMatrix(ddineq);
1052
1053	/Write the normals into class facet/
1054	facet *fAct; //pointer to active facet
1055	int numNonFlip=0;
1056	for (int kk = 0; kk<ddrows; kk++)
1057	{
1058	int ggT=1;//NOTE Why does (int)mpq_get_d(ddineq->matrix[kk][1]) not work?
1059	intvec *load = new intvec(this->numVars);//intvec to store a single facet normal that will then be stored via setFacetNormal
1060	for (int jj = 1; jj <ddcols; jj++)
1061	{
1062	double foo;
1063	foo = mpq_get_d(ddineq->matrix[kk][jj]);
1064	(*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
1065	ggT = intgcd(ggT,(int&)foo);
1066	}//for (int jj = 1; jj <ddcols; jj++)
1067	if(ggT>1)
1068	{
1069	for(int ll=0;ll<this->numVars;ll++)
1070	(*load)[ll] /= ggT;//make primitive vector
1071	}
1072	/*Quick'n'dirty hack for flippability. Executed only if gcone::hasHomInput==FALSE
1073	* Otherwise every facet intersects the positive orthant
1074	*/
1075	if(gcone::hasHomInput==FALSE)
1076	{
1077	//TODO: No dP needed
1078	bool isFlip=FALSE;
1079	for(int jj = 0; jj<load->length(); jj++)
1080	{
1081	// intvec *ivCanonical = new intvec(load->length());
1082	// (*ivCanonical)[jj]=1;
1083	// if (dotProduct(load,ivCanonical)<0)
1084	// {
1085	// isFlip=TRUE;
1086	// break; //URGHS
1087	// }
1088	// delete ivCanonical;
1089	if((*load)[jj]<0)
1090	{
1091	isFlip=TRUE;
1092	break;
1093	}
1094	}/End of check for flippability/
1095	// if(iv64isStrictlyPositive(load))
1096	// isFlip=TRUE;
1097	if(isFlip==FALSE)
1098	{
1099	this->numFacets++;
1100	numNonFlip++;
1101	if(this->numFacets==1)
1102	{
1103	facet *fRoot = new facet();
1104	this->facetPtr = fRoot;
1105	fAct = fRoot;
1106	}
1107	else
1108	{
1109	fAct->next = new facet();
1110	fAct = fAct->next;
1111	}
1112	fAct->isFlippable=FALSE;
1113	fAct->setFacetNormal(load);
1114	fAct->setUCN(this->getUCN());
1115	#ifdef gfan_DEBUG
1116	printf("Marking facet (");load->show(1,0);printf(") as non flippable\n");
1117	#endif
1118	}
1119	else
1120	{
1121	this->numFacets++;
1122	if(this->numFacets==1)
1123	{
1124	facet *fRoot = new facet();
1125	this->facetPtr = fRoot;
1126	fAct = fRoot;
1127	}
1128	else
1129	{
1130	fAct->next = new facet();
1131	fAct = fAct->next;
1132	}
1133	fAct->isFlippable=TRUE;
1134	fAct->setFacetNormal(load);
1135	fAct->setUCN(this->getUCN());
1136	}
1137	}//hasHomInput==FALSE
1138	else //Every facet is flippable
1139	{ /Now load should be full and we can call setFacetNormal/
1140	this->numFacets++;
1141	if(this->numFacets==1)
1142	{
1143	facet *fRoot = new facet();
1144	this->facetPtr = fRoot;
1145	fAct = fRoot;
1146	}
1147	else
1148	{
1149	fAct->next = new facet();
1150	fAct = fAct->next;
1151	}
1152	fAct->isFlippable=TRUE;
1153	fAct->setFacetNormal(load);
1154	fAct->setUCN(this->getUCN());
1155	}//if (isFlippable==FALSE)
1156	delete load;
1157	}//for (int kk = 0; kk<ddrows; kk++)
1158
1159	//In cases like I=<x-1,y-1> there are only non-flippable facets...
1160	if(numNonFlip==this->numFacets)
1161	{
1162	WerrorS ("Only non-flippable facets. Terminating...\n");
1163	// exit(-1);//Bit harsh maybe...
1164	}
1165
1166	/*
1167	Now we should have a linked list containing the facet normals of those facets that are
1168	-irredundant
1169	-flipable
1170	Adressing is done via *facetPtr
1171	*/
1172	if (compIntPoint==TRUE)
1173	{
1174	intvec *iv = new intvec(this->numVars);
1175	dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1176	int jj=1;
1177	for (int ii=0;ii<=this->numVars;ii++)
1178	{
1179	dd_set_si(posRestr->matrix[ii][jj],1);
1180	jj++;
1181	}
1182	dd_MatrixAppendTo(&ddineq,posRestr);
1183	interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets
1184	this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt
1185	delete iv;
1186	dd_FreeMatrix(posRestr);
1187	}
1188	//Clean up but don't delete the return value!
1189	dd_FreeMatrix(ddineq);
1190	set_free(ddredrows);
1191	set_free(ddlinset);
1192	free(ddnewpos);
1193	#ifdef gfanp
1194	gettimeofday(&end, 0);
1195	time_getConeNormals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1196	#endif
1197
1198	}//gcone::getConeNormals(ideal I)
1199
1200	/** \brief Compute the (codim-2)-facets of a given cone
1201	* This method is used during noRevS
1202	* Additionally we check whether the codim2-facet normal is strictly positive. Otherwise
1203	* the facet is marked as non-flippable.
1204	*/
1205	void gcone::getCodim2Normals(const gcone &gc)
1206	{
1207	#ifdef gfanp
1208	timeval start, end;
1209	gettimeofday(&start, 0);
1210	#endif
1211	//this->facetPtr->codim2Ptr = new facet(2); //instantiate a (codim-2)-facet
1212	facet *fAct;
1213	fAct = this->facetPtr;
1214	facet *codim2Act;
1215	//codim2Act = this->facetPtr->codim2Ptr;
1216	dd_MatrixPtr ddineq;//,P,ddakt;
1217	dd_ErrorType err;
1218	//ddineq = facets2Matrix(gc); //get a matrix representation of the cone
1219	ddineq = dd_CopyMatrix(gc.ddFacets);
1220	/Now set appropriate linearity/
1221	for (int ii=0; ii<this->numFacets; ii++)
1222	{
1223	dd_rowset impl_linset, redset;
1224	dd_rowindex newpos;
1225	dd_MatrixPtr ddakt;
1226	ddakt = dd_CopyMatrix(ddineq);
1227	// ddakt->representation=dd_Inequality; //Not using this makes it faster. But why does the quick check below still work?
1228	// ddakt->representation=dd_Generator;
1229	set_addelem(ddakt->linset,ii+1);/Now set appropriate linearity/
1230	#ifdef gfanp
1231	timeval t_ddMC_start, t_ddMC_end;
1232	gettimeofday(&t_ddMC_start,0);
1233	#endif
1234	//dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
1235	dd_PolyhedraPtr ddpolyh;
1236	ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
1237	// ddpolyh=dd_DDMatrix2Poly2(ddakt, dd_MaxCutoff, &err);
1238	dd_MatrixPtr P;
1239	P=dd_CopyGenerators(ddpolyh);
1240	dd_FreePolyhedra(ddpolyh);
1241	//TODO Call for one cone , normalize - check equalities - plus lineality -done
1242	#ifdef gfanp
1243	gettimeofday(&t_ddMC_end,0);
1244	t_ddMC += (t_ddMC_end.tv_sec - t_ddMC_start.tv_sec + 1e-6*(t_ddMC_end.tv_usec - t_ddMC_start.tv_usec));
1245	#endif
1246	/* We loop through each row of P normalize it by making all
1247	* entries integer ones and add the resulting vector to the
1248	* int matrix facet::codim2Facets */
1249	for (int jj=1;jj<=/ddakt/P->rowsize;jj++)
1250	{
1251	fAct->numCodim2Facets++;
1252	if(fAct->numCodim2Facets==1)
1253	{
1254	fAct->codim2Ptr = new facet(2);
1255	codim2Act = fAct->codim2Ptr;
1256	}
1257	else
1258	{
1259	codim2Act->next = new facet(2);
1260	codim2Act = codim2Act->next;
1261	}
1262	intvec *n = new intvec(this->numVars);
1263	#ifdef gfanp
1264	timeval t_mI_start, t_mI_end;
1265	gettimeofday(&t_mI_start,0);
1266	#endif
1267	makeInt(P,jj,*n);
1268	/*for(int kk=0;kk<this->numVars;kk++)
1269	{
1270	int foo;
1271	foo = (int)mpq_get_d(ddakt->matrix[ii][kk+1]);
1272	(*n)[kk]=foo;
1273	}*/
1274	#ifdef gfanp
1275	gettimeofday(&t_mI_end,0);
1276	t_mI += (t_mI_end.tv_sec - t_mI_start.tv_sec + 1e-6*(t_mI_end.tv_usec - t_mI_start.tv_usec));
1277	#endif
1278	codim2Act->setFacetNormal(n);
1279	delete n;
1280	}
1281	/*We check whether the facet spanned by the codim-2 facets
1282	* intersects with the positive orthant. Otherwise we define this
1283	* facet to be non-flippable. Works since we set the appropriate
1284	* linearity for ddakt above.
1285	*/
1286	//TODO It might be faster to compute jus the implied equations instead of a relative interior point
1287	// intvec *iv_intPoint = new intvec(this->numVars);
1288	// dd_MatrixPtr shiftMatrix;
1289	// dd_MatrixPtr intPointMatrix;
1290	// shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
1291	// for(int kk=0;kk<this->numVars;kk++)
1292	// {
1293	// dd_set_si(shiftMatrix->matrix[kk][0],1);
1294	// dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
1295	// }
1296	// intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
1297	// #ifdef gfanp
1298	// timeval t_iP_start, t_iP_end;
1299	// gettimeofday(&t_iP_start, 0);
1300	// #endif
1301	// interiorPoint(intPointMatrix,*iv_intPoint);
1302	// // dd_rowset impl_linste,lbasis;
1303	// // dd_LPSolutionPtr lps=NULL;
1304	// // dd_ErrorType err;
1305	// // dd_FindRelativeInterior(intPointMatrix, &impl_linset, &lbasis, &lps, &err);
1306	// #ifdef gfanp
1307	// gettimeofday(&t_iP_end, 0);
1308	// t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
1309	// #endif
1310	// for(int ll=0;ll<this->numVars;ll++)
1311	// {
1312	// if( (*iv_intPoint)[ll] < 0 )
1313	// {
1314	// fAct->isFlippable=FALSE;
1315	// break;
1316	// }
1317	// }
1318	/End of check/
1319	/This test should be way less time consuming/
1320	#ifdef gfanp
1321	timeval t_iP_start, t_iP_end;
1322	gettimeofday(&t_iP_start, 0);
1323	#endif
1324	bool containsStrictlyPosRay=TRUE;
1325	for(int ii=0;ii<ddakt->rowsize;ii++)
1326	{
1327	containsStrictlyPosRay=TRUE;
1328	for(int jj=1;jj<this->numVars;jj++)
1329	{
1330	if(ddakt->matrix[ii][jj]<=0)
1331	{
1332	containsStrictlyPosRay=FALSE;
1333	break;
1334	}
1335	}
1336	if(containsStrictlyPosRay==TRUE)
1337	break;
1338	}
1339	if(containsStrictlyPosRay==FALSE)
1340	//TODO Not sufficient. Intersect with pos orthant for pos int
1341	fAct->isFlippable=FALSE;
1342	#ifdef gfanp
1343	gettimeofday(&t_iP_end, 0);
1344	t_iP += (t_iP_end.tv_sec - t_iP_start.tv_sec + 1e-6*(t_iP_end.tv_usec - t_iP_start.tv_usec));
1345	#endif
1346	/**/
1347	fAct = fAct->next;
1348	dd_FreeMatrix(ddakt);
1349	dd_FreeMatrix(P);
1350	}//for
1351	dd_FreeMatrix(ddineq);
1352	#ifdef gfanp
1353	gettimeofday(&end, 0);
1354	time_getCodim2Normals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1355	#endif
1356	}
1357
1358	/** Really extremal rays this time ;)
1359	* Extremal rays are unique modulo the homogeneity space.
1360	* Therefore we dd_MatrixAppend gc->ddFacets and gcone::dd_LinealitySpace
1361	* into ddineq. Next we compute the extremal rays of the so given subspace.
1362	* Figuring out whether a ray belongs to a given facet(normal) is done by
1363	* checking whether the inner product of the ray with the normal is zero.
1364	* We use ivAdd here which returns a new intvec. Therefore we need to avoid
1365	* a memory leak which would be cause by the line
1366	* iv=ivAdd(iv,b)
1367	* So we keep pointer tmp to iv and delete(tmp), so there should not occur a
1368	* memleak
1369	* TODO normalization
1370	*/
1371	void gcone::getExtremalRays(const gcone &gc)
1372	{
1373	#ifdef gfanp
1374	timeval start, end;
1375	gettimeofday(&start, 0);
1376	timeval poly_start, poly_end;
1377	gettimeofday(&poly_start,0);
1378	#endif
1379	//Add lineality space - dd_LinealitySpace
1380	dd_MatrixPtr ddineq;
1381	dd_ErrorType err;
1382	if(dd_LinealitySpace->rowsize>0)//The linspace might be 0
1383	ddineq = dd_AppendMatrix(gc.ddFacets,gcone::dd_LinealitySpace);
1384	else
1385	ddineq = dd_CopyMatrix(gc.ddFacets);
1386	/* In case the input is non-homogeneous we add constrains for the positive orthant.
1387	* This is justified by the fact that for non-homog ideals we only consider the
1388	* restricted fan. This way we can be sure to find strictly positive interior points.
1389	* This in turn makes life easy when checking for flippability!
1390	* Drawback: Makes the LP larger so probably slows down computations a wee bit.
1391	*/
1392	dd_MatrixPtr ddPosRestr;
1393	if(hasHomInput==FALSE)
1394	{
1395	dd_MatrixPtr tmp;
1396	ddPosRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1397	for(int ii=0;ii<this->numVars;ii++)
1398	dd_set_si(ddPosRestr->matrix[ii][ii+1],1);
1399	dd_MatrixAppendTo(&ddineq,ddPosRestr);
1400	assert(ddineq);
1401	dd_FreeMatrix(ddPosRestr);
1402	}
1403	dd_PolyhedraPtr ddPolyh;
1404	ddPolyh = dd_DDMatrix2Poly(ddineq, &err);
1405	dd_MatrixPtr P;
1406	P=dd_CopyGenerators(ddPolyh);
1407	dd_FreePolyhedra(ddPolyh);
1408	dd_FreeMatrix(ddineq);
1409	#ifdef gfanp
1410	gettimeofday(&poly_end,0);
1411	t_ddPolyh += (poly_end.tv_sec - poly_start.tv_sec + 1e-6*(poly_end.tv_usec - poly_start.tv_usec));
1412	#endif
1413	/* Compute interior point on the fly*/
1414	intvec *ivIntPointOfCone = new intvec(this->numVars);
1415	mpq_t *colSum = new mpq_t[this->numVars];
1416	int denom[this->numVars];//denominators of colSum
1417	//NOTE TODO need to gcd of rows and factor out! -> makeInt
1418	for(int jj=0;jj<this->numVars;jj++)
1419	{
1420	mpq_init(colSum[jj]);
1421	for(int ii=0;ii<P->rowsize;ii++)
1422	{
1423	mpq_t tmp; mpq_init(tmp);
1424	mpq_t sum; mpq_init(sum);
1425	mpq_set(sum,colSum[jj]);
1426	mpq_add(tmp,sum,P->matrix[ii][jj+1]);
1427	mpq_set(colSum[jj],tmp);
1428	mpq_clear(tmp);
1429	mpq_clear(sum);
1430	}
1431	mpz_t den; mpz_init(den);
1432	mpq_get_den(den,colSum[jj]);
1433	denom[jj]=(int)mpz_get_si(den);
1434	mpz_clear(den);
1435	}
1436	//Now compute lcm of denominators of colSum[jj]
1437	//NOTE gcd as well and normalise instantly?
1438	mpz_t kgV; mpz_init(kgV);
1439	mpz_set_str(kgV,"1",10);
1440	mpz_t den; mpz_init(den);
1441	mpz_t tmp; mpz_init(tmp);
1442	mpq_get_den(tmp,colSum[0]);
1443	for (int ii=0;ii<(this->numVars)-1;ii++)
1444	{
1445	mpq_get_den(den,colSum[ii+1]);
1446	mpz_lcm(kgV,tmp,den);
1447	mpz_set(tmp, kgV);
1448	}
1449	mpq_t qkgV;
1450	mpq_init(qkgV);
1451	mpq_set_z(qkgV,kgV);
1452	mpz_clear(kgV);
1453	mpz_clear(den);
1454	mpz_clear(tmp);
1455	int ggT=(*ivIntPointOfCone)[0];
1456	for (int ii=0;ii<(this->numVars);ii++)
1457	{
1458	mpq_t product;
1459	mpq_init(product);
1460	mpq_mul(product,qkgV,colSum[ii]);
1461	(*ivIntPointOfCone)[ii]=(int)mpz_get_d(mpq_numref(product));
1462	mpq_clear(product);
1463	//Compute intgcd
1464	ggT=intgcd(ggT,(*ivIntPointOfCone)[ii]);
1465	}
1466	//Divide out a common gcd > 1
1467	if(ggT>1)
1468	{
1469	for(int ii=0;ii<this->numVars;ii++)
1470	(*ivIntPointOfCone)[ii] /= ggT;
1471	}
1472	mpq_clear(qkgV);
1473	delete [] colSum;
1474	/For homogeneous input (like Det3,3,5) the int points may be negative. So add a suitable multiple of (1,_,1)/
1475	if(hasHomInput==TRUE && iv64isStrictlyPositive(ivIntPointOfCone)==FALSE)
1476	{
1477	intvec *ivOne = new intvec(this->numVars);
1478	int maxNegEntry=0;
1479	for(int ii=0;ii<this->numVars;ii++)
1480	{
1481	// (*ivOne)[ii]=1;
1482	if ((ivIntPointOfCone)[ii]<maxNegEntry) maxNegEntry=(ivIntPointOfCone)[ii];
1483	}
1484	maxNegEntry *= -1;
1485	maxNegEntry++;//To be on the safe side
1486	for(int ii=0;ii<this->numVars;ii++)
1487	(*ivOne)[ii]=maxNegEntry;
1488	intvec *tmp=ivIntPointOfCone;
1489	ivIntPointOfCone=ivAdd(ivIntPointOfCone,ivOne);
1490	delete(tmp);
1491	// while( !iv64isStrictlyPositive(ivIntPointOfCone) )
1492	// {
1493	// intvec *tmp = ivIntPointOfCone;
1494	// for(int jj=0;jj<this->numVars;jj++)
1495	// (ivOne)[jj] = (ivOne)[jj] << 1; //times 2
1496	// ivIntPointOfCone = ivAdd(ivIntPointOfCone,ivOne);
1497	// delete tmp;
1498	// }
1499	delete ivOne;
1500	int ggT=(*ivIntPointOfCone)[0];
1501	for(int ii=0;ii<this->numVars;ii++)
1502	ggT=intgcd( ggT, (*ivIntPointOfCone)[ii]);
1503	if(ggT>1)
1504	{
1505	for(int jj=0;jj<this->numVars;jj++)
1506	(*ivIntPointOfCone)[jj] /= ggT;
1507	}
1508	}
1509	assert(iv64isStrictlyPositive(ivIntPointOfCone));
1510
1511	this->setIntPoint(ivIntPointOfCone);
1512	delete(ivIntPointOfCone);
1513	/* end of interior point computation*/
1514
1515	//Loop through the rows of P and check whether fNormal*row[i]=0 => row[i] belongs to fNormal
1516	int rows=P->rowsize;
1517	facet *fAct=gc.facetPtr;
1518	//Construct an array to hold the extremal rays of the cone
1519	this->gcRays = (intvec*)omAlloc0(sizeof(intvec)*P->rowsize);
1520	for(int ii=0;ii<P->rowsize;ii++)
1521	{
1522	intvec *rowvec = new intvec(this->numVars);
1523	makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve
1524	this->gcRays[ii] = ivCopy(rowvec);
1525	delete rowvec;
1526	}
1527	this->numRays=P->rowsize;
1528	//Check which rays belong to which facet
1529	while(fAct!=NULL)
1530	{
1531	const intvec *fNormal;// = new intvec(this->numVars);
1532	fNormal = fAct->getRef2FacetNormal();//->getFacetNormal();
1533	intvec *ivIntPointOfFacet = new intvec(this->numVars);
1534	for(int ii=0;ii<rows;ii++)
1535	{
1536	if(dotProduct(*fNormal,this->gcRays[ii])==0)
1537	{
1538	intvec *tmp = ivIntPointOfFacet;//Prevent memleak
1539	fAct->numCodim2Facets++;
1540	facet *codim2Act;
1541	if(fAct->numCodim2Facets==1)
1542	{
1543	fAct->codim2Ptr = new facet(2);
1544	codim2Act = fAct->codim2Ptr;
1545	}
1546	else
1547	{
1548	codim2Act->next = new facet(2);
1549	codim2Act = codim2Act->next;
1550	}
1551	//codim2Act->setFacetNormal(rowvec);
1552	//Rather just let codim2Act point to the corresponding intvec of gcRays
1553	codim2Act->fNormal=this->gcRays[ii];
1554	fAct->numRays++;
1555	//Memleak avoided via tmp
1556	ivIntPointOfFacet=ivAdd(ivIntPointOfFacet,this->gcRays[ii]);
1557	//Now tmp still points to the OLD address of ivIntPointOfFacet
1558	delete(tmp);
1559
1560	}
1561	}//For non-homog input ivIntPointOfFacet should already be >0 here
1562	// if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
1563	//if we have no strictly pos ray but the input is homogeneous
1564	//then add a suitable multiple of (1,...,1)
1565	if( !iv64isStrictlyPositive(ivIntPointOfFacet) && hasHomInput==TRUE)
1566	{
1567	intvec *ivOne = new intvec(this->numVars);
1568	for(int ii=0;ii<this->numVars;ii++)
1569	(*ivOne)[ii]=1;
1570	while( !iv64isStrictlyPositive(ivIntPointOfFacet) )
1571	{
1572	intvec *tmp = ivIntPointOfFacet;
1573	for(int jj=0;jj<this->numVars;jj++)
1574	{
1575	(ivOne)[jj] = (ivOne)[jj] << 1; //times 2
1576	}
1577	ivIntPointOfFacet = ivAdd(ivIntPointOfFacet/diff/,ivOne);
1578	delete tmp;
1579	}
1580	delete ivOne;
1581	}
1582	int ggT=(*ivIntPointOfFacet)[0];
1583	for(int ii=0;ii<this->numVars;ii++)
1584	ggT=intgcd(ggT,(*ivIntPointOfFacet)[ii]);
1585	if(ggT>1)
1586	{
1587	for(int ii=0;ii<this->numVars;ii++)
1588	(*ivIntPointOfFacet)[ii] /= ggT;
1589	}
1590	fAct->setInteriorPoint(ivIntPointOfFacet);
1591
1592	delete(ivIntPointOfFacet);
1593	//Now (if we have at least 3 variables) do a bubblesort on the rays
1594	/*if(this->numVars>2)
1595	{
1596	facet *A[fAct->numRays-1];
1597	facet *f2Act=fAct->codim2Ptr;
1598	for(unsigned ii=0;ii<fAct->numRays;ii++)
1599	{
1600	A[ii]=f2Act;
1601	f2Act=f2Act->next;
1602	}
1603	bool exchanged=FALSE;
1604	unsigned n=fAct->numRays-1;
1605	do
1606	{
1607	exchanged=FALSE;//n=fAct->numRays-1;
1608	for(unsigned ii=0;ii<=n-1;ii++)
1609	{
1610	if((A[ii]->fNormal)->compare((A[ii+1]->fNormal))==1)
1611	{
1612	//Swap rays
1613	cout << "Swapping ";
1614	A[ii]->fNormal->show(1,0); cout << " with "; A[ii+1]->fNormal->show(1,0); cout << endl;
1615	A[ii]->next=A[ii+1]->next;
1616	if(ii>0)
1617	A[ii-1]->next=A[ii+1];
1618	A[ii+1]->next=A[ii];
1619	if(ii==0)
1620	fAct->codim2Ptr=A[ii+1];
1621	//end swap
1622	facet *tmp=A[ii];//swap in list
1623	A[ii+1]=A[ii];
1624	A[ii]=tmp;
1625	// tmp=NULL;
1626	}
1627	}
1628	n--;
1629	}while(exchanged==TRUE && n>=0);
1630	}*///if pVariables>2
1631	// delete fNormal;
1632	fAct = fAct->next;
1633	}//end of facet checking
1634	dd_FreeMatrix(P);
1635	//Now all extremal rays should be set w.r.t their respective fNormal
1636	//TODO Not sufficient -> vol2 II/125&127
1637	//NOTE Sufficient according to cddlibs doc. These ARE rays
1638	//What the hell... let's just take interior points
1639	if(gcone::hasHomInput==FALSE)
1640	{
1641	fAct=gc.facetPtr;
1642	while(fAct!=NULL)
1643	{
1644	// bool containsStrictlyPosRay=FALSE;
1645	// facet *codim2Act;
1646	// codim2Act = fAct->codim2Ptr;
1647	// while(codim2Act!=NULL)
1648	// {
1649	// intvec *rayvec;
1650	// rayvec = codim2Act->getFacetNormal();//Mind this is no normal but a ray!
1651	// //int negCtr=0;
1652	// if(iv64isStrictlyPositive(rayvec))
1653	// {
1654	// containsStrictlyPosRay=TRUE;
1655	// delete(rayvec);
1656	// break;
1657	// }
1658	// delete(rayvec);
1659	// codim2Act = codim2Act->next;
1660	// }
1661	// if(containsStrictlyPosRay==FALSE)
1662	// fAct->isFlippable=FALSE;
1663	if(!iv64isStrictlyPositive(fAct->interiorPoint))
1664	fAct->isFlippable=FALSE;
1665	fAct = fAct->next;
1666	}
1667	}//hasHomInput?
1668	#ifdef gfanp
1669	gettimeofday(&end, 0);
1670	t_getExtremalRays += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
1671	#endif
1672	}
1673
1674	inline bool gcone::iv64isStrictlyPositive(const intvec * iv64)
1675	{
1676	bool res=TRUE;
1677	for(int ii=0;ii<iv64->length();ii++)
1678	{
1679	if((*iv64)[ii]<=0)
1680	{
1681	res=FALSE;
1682	break;
1683	}
1684	}
1685	return res;
1686	}
1687
1688	/** \brief Compute the Groebner Basis on the other side of a shared facet
1689	*
1690	* Implements algorithm 4.3.2 from Anders' thesis.
1691	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
1692	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
1693	* is parallel to \f$ leadexp(g) \f$
1694	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
1695	* Other possibilities include computing the rank of the matrix consisting of the vectors in question and
1696	* computing an interior point of the facet and taking all terms having the same weight with respect
1697	* to this interior point.
1698	*\param ideal, facet
1699	* Input: a marked,reduced Groebner basis and a facet
1700	*/
1701	inline void gcone::flip(ideal gb, facet *f) //Compute "the other side"
1702	{
1703	#ifdef gfanp
1704	timeval start, end;
1705	gettimeofday(&start, 0);
1706	#endif
1707	intvec *fNormal;// = new intvec(this->numVars); //facet normal, check for parallelity
1708	fNormal = f->getFacetNormal(); //read this->fNormal;
1709	#ifdef gfan_DEBUG
1710	// std::cout << "running gcone::flip" << std::endl;
1711	printf("flipping UCN %i over facet",this->getUCN());
1712	fNormal->show(1,0);
1713	printf(") with UCN %i\n",f->getUCN() );
1714	#endif
1715	if(this->getUCN() != f->getUCN())
1716	{
1717	WerrorS("Uh oh... Trying to flip over facet with incompatible UCN");
1718	exit(-1);
1719	}
1720	/1st step: Compute the initial ideal/
1721	/poly initialFormElement[IDELEMS(gb)];/ //array of #polys in GB to store initial form
1722	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
1723
1724	computeInv(gb,initialForm,*fNormal);
1725
1726	#ifdef gfan_DEBUG
1727	/* cout << "Initial ideal is: " << endl;
1728	idShow(initialForm);
1729	//f->printFlipGB();*/
1730	// cout << "===" << endl;
1731	#endif
1732	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
1733	/*Substep 2.1
1734	compute $G_{-\alpha}(in_v(I))
1735	see journal p. 66
1736	NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
1737	srcRing already has a weighted ordering
1738	*/
1739	ring srcRing=currRing;
1740	ring tmpRing;
1741
1742	if( (srcRing->order[0]!=ringorder_a))
1743	{
1744	intvec *iv;// = new intvec(this->numVars);
1745	iv = ivNeg(fNormal);//ivNeg uses ivCopy -> new
1746	// tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
1747	tmpRing=rCopyAndAddWeight(srcRing,iv);
1748	delete iv;
1749	}
1750	else
1751	{
1752	tmpRing=rCopy0(srcRing);
1753	int length=fNormal->length();
1754	int A=(int )omAlloc0(length*sizeof(int));
1755	for(int jj=0;jj<length;jj++)
1756	{
1757	A[jj]=-(*fNormal)[jj];
1758	}
1759	omFree(tmpRing->wvhdl[0]);
1760	tmpRing->wvhdl[0]=(int*)A;
1761	tmpRing->block1[0]=length;
1762	rComplete(tmpRing);
1763	//omFree(A);
1764	}
1765	delete fNormal;
1766	rChangeCurrRing(tmpRing);
1767
1768	ideal ina;
1769	ina=idrCopyR(initialForm,srcRing);
1770	idDelete(&initialForm);
1771	ideal H;
1772	// H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
1773	#ifdef gfanp
1774	timeval t_kStd_start, t_kStd_end;
1775	gettimeofday(&t_kStd_start,0);
1776	#endif
1777	if(gcone::hasHomInput==TRUE)
1778	H=kStd(ina,NULL,isHomog,NULL/,gcone::hilbertFunction/);
1779	else
1780	H=kStd(ina,NULL,isNotHomog,NULL); //This is \mathcal(G)_{>_-\alpha}(in_v(I))
1781	#ifdef gfanp
1782	gettimeofday(&t_kStd_end, 0);
1783	t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
1784	#endif
1785	idSkipZeroes(H);
1786	idDelete(&ina);
1787
1788	/*Substep 2.2
1789	do the lifting and mark according to H
1790	*/
1791	rChangeCurrRing(srcRing);
1792	ideal srcRing_H;
1793	ideal srcRing_HH;
1794	srcRing_H=idrCopyR(H,tmpRing);
1795	//H is needed further below, so don't idDelete here
1796	srcRing_HH=ffG(srcRing_H,this->gcBasis);
1797	idDelete(&srcRing_H);
1798
1799	/*Substep 2.2.1
1800	* Mark according to G_-\alpha
1801	* Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
1802	* we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
1803	* represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
1804	* Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
1805	* compute the difference accordingly
1806	*/
1807	#ifdef gfanp
1808	timeval t_markings_start, t_markings_end;
1809	gettimeofday(&t_markings_start, 0);
1810	#endif
1811	bool markingsAreCorrect=FALSE;
1812	dd_MatrixPtr intPointMatrix;
1813	int iPMatrixRows=0;
1814	dd_rowrange aktrow=0;
1815	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1816	{
1817	poly aktpoly=(poly)srcRing_HH->m[ii];//This is a pointer, so don't pDelete
1818	iPMatrixRows = iPMatrixRows+pLength(aktpoly);
1819	}
1820	/* additionally one row for the standard-simplex and another for a row that becomes 0 during
1821	* construction of the differences
1822	*/
1823	intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
1824	intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
1825
1826	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1827	{
1828	markingsAreCorrect=FALSE; //crucial to initialise here
1829	poly aktpoly=srcRing_HH->m[ii]; //Only a pointer, so don't pDelete
1830	/Comparison of leading monomials is done via exponent vectors/
1831	for (int jj=0;jj<IDELEMS(H);jj++)
1832	{
1833	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
1834	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
1835	pGetExpV(aktpoly,src_ExpV);
1836	rChangeCurrRing(tmpRing); //this ring change is crucial!
1837	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
1838	rChangeCurrRing(srcRing);
1839	bool expVAreEqual=TRUE;
1840	for (int kk=1;kk<=this->numVars;kk++)
1841	{
1842	#ifdef gfan_DEBUG
1843	// cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
1844	#endif
1845	if (src_ExpV[kk]!=dst_ExpV[kk])
1846	{
1847	expVAreEqual=FALSE;
1848	}
1849	}
1850	if (expVAreEqual==TRUE)
1851	{
1852	markingsAreCorrect=TRUE; //everything is fine
1853	#ifdef gfan_DEBUG
1854	// cout << "correct markings" << endl;
1855	#endif
1856	}//if (pHead(aktpoly)==pHead(H->m[jj])
1857	omFree(src_ExpV);
1858	omFree(dst_ExpV);
1859	}//for (int jj=0;jj<IDELEMS(H);jj++)
1860
1861	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
1862	if (markingsAreCorrect==TRUE)
1863	{
1864	pGetExpV(aktpoly,leadExpV);
1865	}
1866	else
1867	{
1868	rChangeCurrRing(tmpRing);
1869	pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
1870	rChangeCurrRing(srcRing);
1871	}
1872	/compute differences of the expvects/
1873	while (pNext(aktpoly)!=NULL)
1874	{
1875	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1876	/*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
1877	is not omitted when computing the differences*/
1878	if(markingsAreCorrect==TRUE)
1879	{
1880	aktpoly=pNext(aktpoly);
1881	pGetExpV(aktpoly,v);
1882	}
1883	else
1884	{
1885	pGetExpV(pHead(aktpoly),v);
1886	markingsAreCorrect=TRUE;
1887	}
1888	int ctr=0;
1889	for (int jj=0;jj<this->numVars;jj++)
1890	{
1891	/Store into ddMatrix/
1892	if(leadExpV[jj+1]-v[jj+1])
1893	ctr++;
1894	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
1895	}
1896	/It ought to be more sensible to avoid 0-rows in the first place/
1897	// if(ctr==this->numVars)//We have a 0-row
1898	// dd_MatrixRowRemove(&intPointMatrix,aktrow);
1899	// else
1900	aktrow +=1;
1901	omFree(v);
1902	}
1903	omFree(leadExpV);
1904	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1905	#ifdef gfanp
1906	gettimeofday(&t_markings_end, 0);
1907	t_markings += (t_markings_end.tv_sec - t_markings_start.tv_sec + 1e-6*(t_markings_end.tv_usec - t_markings_start.tv_usec));
1908	#endif
1909	/Now it is safe to idDelete(H)/
1910	idDelete(&H);
1911	/*Preprocessing goes here since otherwise we would delete the constraint
1912	* for the standard simplex.
1913	*/
1914	preprocessInequalities(intPointMatrix);
1915	/Now we add the constraint for the standard simplex/
1916	// dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
1917	// for (int jj=1;jj<=this->numVars;jj++)
1918	// {
1919	// dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
1920	// }
1921	//Let's make sure we compute interior points from the positive orthant
1922	// dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1923	//
1924	// int jj=1;
1925	// for (int ii=0;ii<this->numVars;ii++)
1926	// {
1927	// dd_set_si(posRestr->matrix[ii][jj],1);
1928	// jj++;
1929	// }
1930	/*We create a matrix containing the standard simplex
1931	* and constraints to assure a strictly positive point
1932	* is computed */
1933	dd_MatrixPtr posRestr = dd_CreateMatrix(this->numVars+1, this->numVars+1);
1934	for(int ii=0;ii<posRestr->rowsize;ii++)
1935	{
1936	if(ii==0)
1937	{
1938	dd_set_si(posRestr->matrix[ii][0],-1);
1939	for(int jj=1;jj<=this->numVars;jj++)
1940	dd_set_si(posRestr->matrix[ii][jj],1);
1941	}
1942	else
1943	{
1944	/** Set all variables to \geq 1/10. YMMV but this choice is pretty equal*/
1945	dd_set_si2(posRestr->matrix[ii][0],-1,2);
1946	dd_set_si(posRestr->matrix[ii][ii],1);
1947	}
1948	}
1949	dd_MatrixAppendTo(&intPointMatrix,posRestr);
1950	dd_FreeMatrix(posRestr);
1951
1952	intvec *iv_weight = new intvec(this->numVars);
1953	#ifdef gfanp
1954	timeval t_dd_start, t_dd_end;
1955	gettimeofday(&t_dd_start, 0);
1956	#endif
1957	dd_ErrorType err;
1958	dd_rowset implLin, redrows;
1959	dd_rowindex newpos;
1960
1961	//NOTE Here we should remove interiorPoint and instead
1962	// create and ordering like (a(omega),a(fNormal),dp)
1963	// if(this->ivIntPt==NULL)
1964	interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point
1965	// else
1966	// iv_weight=this->getIntPoint();
1967	dd_FreeMatrix(intPointMatrix);
1968	/Crude attempt for interior point /
1969	/*dd_PolyhedraPtr ddpolyh;
1970	dd_ErrorType err;
1971	dd_rowset impl_linset,redset;
1972	dd_rowindex newpos;
1973	dd_MatrixCanonicalize(&intPointMatrix, &impl_linset, &redset, &newpos, &err);
1974	ddpolyh=dd_DDMatrix2Poly(intPointMatrix, &err);
1975	dd_MatrixPtr P;
1976	P=dd_CopyGenerators(ddpolyh);
1977	dd_FreePolyhedra(ddpolyh);
1978	for(int ii=0;ii<P->rowsize;ii++)
1979	{
1980	intvec *iv_row=new intvec(this->numVars);
1981	makeInt(P,ii+1,*iv_row);
1982	iv_weight =ivAdd(iv_weight, iv_row);
1983	delete iv_row;
1984	}
1985	dd_FreeMatrix(P);
1986	dd_FreeMatrix(intPointMatrix);*/
1987	#ifdef gfanp
1988	gettimeofday(&t_dd_end, 0);
1989	t_dd += (t_dd_end.tv_sec - t_dd_start.tv_sec + 1e-6*(t_dd_end.tv_usec - t_dd_start.tv_usec));
1990	#endif
1991
1992	/*Step 3
1993	* turn the minimal basis into a reduced one */
1994	// NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
1995	// Thus:
1996	//ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
1997	ring dstRing=rCopy0(tmpRing);
1998	int length=iv_weight->length();
1999	int A=(int )omAlloc0(length*sizeof(int));
2000	for(int jj=0;jj<length;jj++)
2001	{
2002	A[jj]=(*iv_weight)[jj];
2003	}
2004	dstRing->wvhdl[0]=(int*)A;
2005	rComplete(dstRing);
2006	rChangeCurrRing(dstRing);
2007	rDelete(tmpRing);
2008	delete iv_weight;
2009
2010	ideal dstRing_I;
2011	dstRing_I=idrCopyR(srcRing_HH,srcRing);
2012	idDelete(&srcRing_HH); //Hmm.... causes trouble - no more
2013	//dstRing_I=idrCopyR(inputIdeal,srcRing);
2014	BITSET save=test;
2015	test\|=Sy_bit(OPT_REDSB);
2016	test\|=Sy_bit(OPT_REDTAIL);
2017	#ifdef gfan_DEBUG
2018	// test\|=Sy_bit(6); //OPT_DEBUG
2019	#endif
2020	ideal tmpI;
2021	//NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
2022	//tmpI = idrCopyR(this->inputIdeal,this->baseRing);
2023	tmpI = dstRing_I;
2024	#ifdef gfanp
2025	gettimeofday(&t_kStd_start,0);
2026	#endif
2027	if(gcone::hasHomInput==TRUE)
2028	dstRing_I=kStd(tmpI,NULL,isHomog,NULL/,gcone::hilbertFunction/);
2029	else
2030	dstRing_I=kStd(tmpI,NULL,isNotHomog,NULL);
2031	#ifdef gfanp
2032	gettimeofday(&t_kStd_end, 0);
2033	t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
2034	#endif
2035	idDelete(&tmpI);
2036	idNorm(dstRing_I);
2037	// kInterRed(dstRing_I);
2038	idSkipZeroes(dstRing_I);
2039	test=save;
2040	/End of step 3 - reduction/
2041
2042	f->setFlipGB(dstRing_I);//store the flipped GB
2043	// idDelete(&dstRing_I);
2044	f->flipRing=rCopy(dstRing); //store the ring on the other side
2045	#ifdef gfan_DEBUG
2046	printf("Flipped GB is UCN %i:\n",counter+1);
2047	idDebugPrint(dstRing_I);
2048	printf("\n");
2049	#endif
2050	idDelete(&dstRing_I);
2051	rChangeCurrRing(srcRing); //return to the ring we started the computation of flipGB in
2052	rDelete(dstRing);
2053	#ifdef gfanp
2054	gettimeofday(&end, 0);
2055	time_flip += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2056	#endif
2057	}//void flip(ideal gb, facet *f)
2058
2059	/** \brief A slightly different approach to flipping
2060	* Here we use the fact that in_v(in_u(I))=in_(u+eps*v)(I). Therefore, we do no longer
2061	* need to compute an interior point and run BBA on the minimal basis but we can rather
2062	* use the ordering (a(omega),a(fNormal),dp)
2063	* The second parameter facet *f must not be const since we need to store f->flipGB
2064	* Problem: Assume we start in a cone with ordering (dp,C). Then \f$ in_\omega(I) \f$
2065	* will be from a ring with (a(),dp,C) and our resulting cone from (a(),a(),dp,C). Hence a way
2066	* must be found to circumvent the sequence of a()'s growing to a ridiculous size.
2067	* Therefore: We use (a(),a(),dp,C) for the computation of the reduced basis. But then we
2068	* do have an interior point of the cone by adding the extremal rays. So we replace
2069	* the latter cone by a cone with (a(sum_of_rays),dp,C).
2070	* Con: It's incredibly ugly
2071	* Pro: No messing around with readConeFromFile()
2072	* Is there a way to construct a vector from \f$ \omega \f$ and the facet normal?
2073	*/
2074	inline void gcone::flip2(const ideal &gb, facet *f)
2075	{
2076	#ifdef gfanp
2077	timeval start, end;
2078	gettimeofday(&start, 0);
2079	#endif
2080	const intvec *fNormal;
2081	fNormal = f->getRef2FacetNormal();/->getFacetNormal();/ //read this->fNormal;
2082	#ifdef gfan_DEBUG
2083	printf("flipping UCN %i over facet(",this->getUCN());
2084	fNormal->show(1,0);
2085	printf(") with UCN %i\n",f->getUCN());
2086	#endif
2087	if(this->getUCN() != f->getUCN())
2088	{ printf("%i vs %i\n",this->getUCN(), f->getUCN() );
2089	WerrorS("Uh oh... Trying to flip over facet with incompatible UCN");
2090	exit(-1);
2091	}
2092	/1st step: Compute the initial ideal/
2093	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
2094	computeInv( gb, initialForm, *fNormal );
2095	ring srcRing=currRing;
2096	ring tmpRing;
2097
2098	const intvec *intPointOfFacet;
2099	intPointOfFacet=f->getInteriorPoint();
2100	//Now we need two blocks of ringorder_a!
2101	//May assume the same situation as in flip() here
2102	if( (srcRing->order[0]!=ringorder_a) && (srcRing->order[1]!=ringorder_a) )
2103	{
2104	intvec *iv = new intvec(this->numVars);//init with 1s, since we do not need a 2nd block here but later
2105	// intvec *iv_foo = new intvec(this->numVars,1);//placeholder
2106	intvec *ivw = ivNeg(fNormal);
2107	tmpRing=rCopyAndAddWeight2(srcRing,ivw/intPointOfFacet/,iv);
2108	delete iv;delete ivw;
2109	// delete iv_foo;
2110	}
2111	else
2112	{
2113	intvec *iv=new intvec(this->numVars);
2114	intvec *ivw=ivNeg(fNormal);
2115	tmpRing=rCopyAndAddWeight2(srcRing,ivw,iv);
2116	delete iv; delete ivw;
2117	/*tmpRing=rCopy0(srcRing);
2118	int length=fNormal->length();
2119	int A1=(int )omAlloc0(length*sizeof(int));
2120	int A2=(int )omAlloc0(length*sizeof(int));
2121	for(int jj=0;jj<length;jj++)
2122	{
2123	A1[jj] = -(*fNormal)[jj];
2124	A2[jj] = 1;//-(*fNormal)[jj];//NOTE Do we need this here? This is only the facet ideal
2125	}
2126	omFree(tmpRing->wvhdl[0]);
2127	if(tmpRing->wvhdl[1]!=NULL)
2128	omFree(tmpRing->wvhdl[1]);
2129	tmpRing->wvhdl[0]=(int*)A1;
2130	tmpRing->block1[0]=length;
2131	tmpRing->wvhdl[1]=(int*)A2;
2132	tmpRing->block1[1]=length;
2133	rComplete(tmpRing);*/
2134	}
2135	// delete fNormal; //NOTE Do not delete when using getRef2FacetNormal();
2136	rChangeCurrRing(tmpRing);
2137	//Now currRing should have (a(),a(),dp,C)
2138	ideal ina;
2139	ina=idrCopyR(initialForm,srcRing);
2140	idDelete(&initialForm);
2141	ideal H;
2142	#ifdef gfanp
2143	timeval t_kStd_start, t_kStd_end;
2144	gettimeofday(&t_kStd_start,0);
2145	#endif
2146	BITSET save=test;
2147	test\|=Sy_bit(OPT_REDSB);
2148	test\|=Sy_bit(OPT_REDTAIL);
2149	// if(gcone::hasHomInput==TRUE)
2150	H=kStd(ina,NULL,testHomog/isHomog/,NULL/,gcone::hilbertFunction/);
2151	// else
2152	// H=kStd(ina,NULL,isNotHomog,NULL); //This is \mathcal(G)_{>_-\alpha}(in_v(I))
2153	test=save;
2154	#ifdef gfanp
2155	gettimeofday(&t_kStd_end, 0);
2156	t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
2157	#endif
2158	idSkipZeroes(H);
2159	idDelete(&ina);
2160
2161	rChangeCurrRing(srcRing);
2162	ideal srcRing_H;
2163	ideal srcRing_HH;
2164	srcRing_H=idrCopyR(H,tmpRing);
2165	//H is needed further below, so don't idDelete here
2166	srcRing_HH=ffG(srcRing_H,this->gcBasis);
2167	idDelete(&srcRing_H);
2168	//Now BBA(srcRing_HH) with (a(),a(),dp)
2169	/* Evil modification of currRing */
2170	ring dstRing=rCopy0(tmpRing);
2171	int length=this->numVars;
2172	int A1=(int )omAlloc0(length*sizeof(int));
2173	int A2=(int )omAlloc0(length*sizeof(int));
2174	const intvec *ivw=f->getRef2FacetNormal();
2175	for(int jj=0;jj<length;jj++)
2176	{
2177	A1[jj] = (*intPointOfFacet)[jj];
2178	A2[jj] = -(ivw)[jj];//TODO Or minus (ivw)[ii] ??? NOTE minus
2179	}
2180	omFree(dstRing->wvhdl[0]);
2181	if(dstRing->wvhdl[1]!=NULL)
2182	omFree(dstRing->wvhdl[1]);
2183	dstRing->wvhdl[0]=(int*)A1;
2184	dstRing->block1[0]=length;
2185	dstRing->wvhdl[1]=(int*)A2;
2186	dstRing->block1[1]=length;
2187	rComplete(dstRing);
2188	rChangeCurrRing(dstRing);
2189	ideal dstRing_I;
2190	dstRing_I=idrCopyR(srcRing_HH,srcRing);
2191	idDelete(&srcRing_HH); //Hmm.... causes trouble - no more
2192	save=test;
2193	test\|=Sy_bit(OPT_REDSB);
2194	test\|=Sy_bit(OPT_REDTAIL);
2195	ideal tmpI;
2196	tmpI = dstRing_I;
2197	#ifdef gfanp
2198	// timeval t_kStd_start, t_kStd_end;
2199	gettimeofday(&t_kStd_start,0);
2200	#endif
2201	// if(gcone::hasHomInput==TRUE)
2202	// dstRing_I=kStd(tmpI,NULL,isHomog,NULL/,gcone::hilbertFunction/);
2203	// else
2204	dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
2205	#ifdef gfanp
2206	gettimeofday(&t_kStd_end, 0);
2207	t_kStd += (t_kStd_end.tv_sec - t_kStd_start.tv_sec + 1e-6*(t_kStd_end.tv_usec - t_kStd_start.tv_usec));
2208	#endif
2209	idDelete(&tmpI);
2210	idNorm(dstRing_I);
2211	idSkipZeroes(dstRing_I);
2212	test=save;
2213	/End of step 3 - reduction/
2214
2215	f->setFlipGB(dstRing_I);
2216	f->flipRing=rCopy(dstRing);
2217	rDelete(tmpRing);
2218	rDelete(dstRing);
2219	//Now we should have dstRing with (a(),a(),dp,C)
2220	//This must be replaced with (a(),dp,C) BEFORE gcTmp is actually added to the list
2221	//of cones in noRevS
2222	rChangeCurrRing(srcRing);
2223	#ifdef gfanp
2224	gettimeofday(&end, 0);
2225	time_flip2 += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2226	#endif
2227	}//flip2
2228
2229	/** \brief Compute initial ideal
2230	* Compute the initial ideal in_v(G) wrt a (possible) facet normal
2231	* used in gcone::getFacetNormal in order to preprocess possible facet normals
2232	* and in gcone::flip for obvious reasons.
2233	*/
2234	inline void gcone::computeInv(const ideal &gb, ideal &initialForm, const intvec &fNormal)
2235	{
2236	#ifdef gfanp
2237	timeval start, end;
2238	gettimeofday(&start, 0);
2239	#endif
2240	for (int ii=0;ii<IDELEMS(gb);ii++)
2241	{
2242	poly initialFormElement;
2243	poly aktpoly = (poly)gb->m[ii];//Ptr, so don't pDelete(aktpoly)
2244	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
2245	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
2246	initialFormElement=pHead(aktpoly);
2247	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
2248	{
2249	intvec *check = new intvec(this->numVars);
2250	aktpoly=pNext(aktpoly); //next term
2251	// pSetm(aktpoly);
2252	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
2253	pGetExpV(aktpoly,v);
2254	/* Convert (int)v into (intvec)check */
2255	for (int jj=0;jj<this->numVars;jj++)
2256	{
2257	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
2258	}
2259	if (isParallel(check,fNormal)) //pass check when
2260	// if(isParallel((const intvec)&check,(const intvec)&fNormal))
2261	// if(fNormal.compare(check)==0)
2262	{
2263	//Found a parallel vector. Add it
2264	initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));
2265	}
2266	omFree(v);
2267	delete check;
2268	}//while
2269	#ifdef gfan_DEBUG
2270	// cout << "Initial Form=";
2271	// pWrite(initialFormElement[ii]);
2272	// cout << "---" << endl;
2273	#endif
2274	/Now initialFormElement must be added to (ideal)initialForm /
2275	initialForm->m[ii]=pCopy(initialFormElement);
2276	pDelete(&initialFormElement);
2277	omFree(leadExpV);
2278	}//for
2279	#ifdef gfanp
2280	gettimeofday(&end, 0);
2281	time_computeInv += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2282	#endif
2283	}
2284
2285	/** \brief Compute the remainder of a polynomial by a given ideal
2286	*
2287	* Compute \f$ f^{\mathcal{G}} \f$
2288	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
2289	* However, since we are only interested in the remainder, there is no need to
2290	* compute the factors \f$ a_i \f$
2291	*/
2292	//NOTE: Should be replaced by kNF or kNF2
2293	//NOTE: Done
2294	//NOTE: removed with r12286
2295
2296	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
2297	*/
2298	//NOTE: use kNF or kNF2 instead of restOfDivision
2299	inline ideal gcone::ffG(const ideal &H, const ideal &G)
2300	{
2301	int size=IDELEMS(H);
2302	ideal res=idInit(size,1);
2303	for (int ii=0;ii<size;ii++)
2304	{
2305	// poly temp1;//=pInit();
2306	// poly temp2;//=pInit();
2307	poly temp3;//=pInit();//polys to temporarily store values for pSub
2308	// res->m[ii]=pCopy(kNF(G, NULL,H->m[ii],0,0));
2309	// temp1=pCopy(H->m[ii]);//TRY
2310	// temp2=pCopy(res->m[ii]);
2311	//NOTE if gfanHeuristic=0 (sic!) this results in dPolyErrors - mon from wrong ring
2312	// temp2=pCopy(kNF(G, NULL,H->m[ii],0,0));//TRY
2313	// temp3=pSub(temp1, temp2);//TRY
2314	temp3=pSub(pCopy(H->m[ii]),pCopy(kNF(G,NULL,H->m[ii],0,0)));//NOTRY
2315	res->m[ii]=pCopy(temp3);
2316	//res->m[ii]=pSub(temp1,temp2); //buggy
2317	//cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
2318	// pDelete(&temp1);//TRY
2319	// pDelete(&temp2);
2320	pDelete(&temp3);
2321	}
2322	return res;
2323	}
2324
2325	/** \brief Preprocessing of inequlities
2326	* Do some preprocessing on the matrix of inequalities
2327	* 1) Replace several constraints on the pos. orthants by just one for each orthant
2328	* 2) Remove duplicates of inequalities
2329	* 3) Remove inequalities that arise as sums of other inequalities
2330	*/
2331	void gcone::preprocessInequalities(dd_MatrixPtr &ddineq)
2332	{
2333	/* int *posRowsArray=NULL;
2334	int num_alloc=0;
2335	int num_elts=0;
2336	int offset=0;*/
2337	//Remove zeroes (and strictly pos rows?)
2338	for(int ii=0;ii<ddineq->rowsize;ii++)
2339	{
2340	intvec *iv = new intvec(this->numVars);
2341	intvec ivNull = new intvec(this->numVars);//Needed for intvec64::compare(intvec)
2342	int posCtr=0;
2343	for(int jj=0;jj<this->numVars;jj++)
2344	{
2345	(*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
2346	if((*iv)[jj]>0)//check for strictly pos rows
2347	posCtr++;
2348	//Behold! This will delete the row for the standard simplex!
2349	}
2350	// if( (iv->compare(0)==0) \|\| (posCtr==iv->length()) )
2351	if( (posCtr==iv->length()) \|\| (iv->compare(ivNull)==0) )
2352	{
2353	dd_MatrixRowRemove(&ddineq,ii+1);
2354	ii--;//Yes. This is on purpose
2355	}
2356	delete iv;
2357	delete ivNull;
2358	}
2359	//Remove duplicates of rows
2360	// posRowsArray=NULL;
2361	// num_alloc=0;
2362	// num_elts=0;
2363	// offset=0;
2364	// int num_newRows = ddineq->rowsize;
2365	// for(int ii=0;ii<ddineq->rowsize-1;ii++)
2366	// for(int ii=0;ii<num_newRows-1;ii++)
2367	// {
2368	// intvec *iv = new intvec(this->numVars);//1st vector to check against
2369	// for(int jj=0;jj<this->numVars;jj++)
2370	// (*iv)[jj]=(int)mpq_get_d(ddineq->matrix[ii][jj+1]);
2371	// for(int jj=ii+1;jj</ddineq->rowsize/num_newRows;jj++)
2372	// {
2373	// intvec *ivCheck = new intvec(this->numVars);//Checked against iv
2374	// for(int kk=0;kk<this->numVars;kk++)
2375	// (*ivCheck)[kk]=(int)mpq_get_d(ddineq->matrix[jj][kk+1]);
2376	// if (iv->compare(ivCheck)==0)
2377	// {
2378	// // cout << "=" << endl;
2379	// // num_alloc++;
2380	// // void tmp=realloc(posRowsArray,(num_allocsizeof(int)));
2381	// // if(!tmp)
2382	// // {
2383	// // WerrorS("Woah dude! Couldn't realloc memory\n");
2384	// // exit(-1);
2385	// // }
2386	// // posRowsArray = (int*)tmp;
2387	// // posRowsArray[num_elts]=jj;
2388	// // num_elts++;
2389	// dd_MatrixRowRemove(&ddineq,jj+1);
2390	// num_newRows = ddineq->rowsize;
2391	// }
2392	// delete ivCheck;
2393	// }
2394	// delete iv;
2395	// }
2396	// for(int ii=0;ii<num_elts;ii++)
2397	// {
2398	// dd_MatrixRowRemove(&ddineq,posRowsArray[ii]+1-offset);
2399	// offset++;
2400	// }
2401	// free(posRowsArray);
2402	//Apply Thm 2.1 of JOTA Vol 53 No 1 April 1987*/
2403	}//preprocessInequalities
2404
2405	/** \brief Compute a Groebner Basis
2406	*
2407	* Computes the Groebner basis and stores the result in gcone::gcBasis
2408	*\param ideal
2409	*\return void
2410	*/
2411	inline void gcone::getGB(const ideal &inputIdeal)
2412	{
2413	BITSET save=test;
2414	test\|=Sy_bit(OPT_REDSB);
2415	test\|=Sy_bit(OPT_REDTAIL);
2416	ideal gb;
2417	gb=kStd(inputIdeal,NULL,testHomog,NULL);
2418	idNorm(gb);
2419	idSkipZeroes(gb);
2420	this->gcBasis=gb; //write the GB into gcBasis
2421	test=save;
2422	}//void getGB
2423
2424	/** \brief Compute the negative of a given intvec
2425	*/
2426	static intvec* ivNeg(const intvec *iv)
2427	{ //Hm, switching to intvec const intvec does no longer work
2428	intvec *res;// = new intvec(iv->length());
2429	res=ivCopy(iv);
2430	res = (int)-1;
2431	return res;
2432	}
2433
2434
2435	/** \brief Compute the dot product of two intvecs
2436	*
2437	*/
2438	static int dotProduct(const intvec &iva, const intvec &ivb)
2439	{
2440	int res=0;
2441	for (int i=0;i<pVariables;i++)
2442	{
2443	// #ifndef NDEBUG
2444	// (const_cast<intvec>(&iva))->show(1,0); (const_cast<intvec>(&ivb))->show(1,0);
2445	// #endif
2446	res = res+(iva[i]*ivb[i]);
2447	}
2448	return res;
2449	}
2450	/** \brief Check whether two intvecs are parallel
2451	*
2452	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
2453	*/
2454	static bool isParallel(const intvec &a,const intvec &b)
2455	{
2456	/*#ifdef gfanp
2457	timeval start, end;
2458	gettimeofday(&start, 0);
2459	#endif*/
2460	int lhs,rhs;
2461	bool res;
2462	lhs=dotProduct(a,b)*dotProduct(a,b);
2463	rhs=dotProduct(a,a)*dotProduct(b,b);
2464	//cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
2465	if (lhs==rhs)
2466	{
2467	res = TRUE;
2468	}
2469	else
2470	{
2471	res = FALSE;
2472	}
2473	// #ifdef gfanp
2474	// gettimeofday(&end, 0);
2475	// t_isParallel += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
2476	// #endif
2477	return res;
2478	}//bool isParallel
2479
2480	/** \brief Compute an interior point of a given cone
2481	* Result will be written into intvec iv.
2482	* Any rational point is automatically converted into an integer.
2483	*/
2484	void gcone::interiorPoint( dd_MatrixPtr &M, intvec &iv) //no const &M here since we want to remove redundant rows
2485	{
2486	dd_LPPtr lp,lpInt;
2487	dd_ErrorType err=dd_NoError;
2488	dd_LPSolverType solver=dd_DualSimplex;
2489	dd_LPSolutionPtr lpSol=NULL;
2490	dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior
2491	dd_rowindex ddnewpos;
2492	dd_NumberType numb;
2493	//M->representation=dd_Inequality;
2494
2495	//NOTE: Make this n-dimensional!
2496	//dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
2497
2498	/*NOTE: Leave the following line commented out!
2499	* Otherwise it will slow down computations a great deal
2500	* */
2501	// dd_MatrixCanonicalizeLinearity(&M, &ddlinset, &ddnewpos, &err);
2502	//if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
2503	dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
2504	int jj=1;
2505	for (int ii=0;ii<this->numVars;ii++)
2506	{
2507	dd_set_si(posRestr->matrix[ii][jj],1);
2508	jj++;
2509	}
2510	dd_MatrixAppendTo(&M,posRestr);
2511	dd_FreeMatrix(posRestr);
2512	lp=dd_Matrix2LP(M, &err);
2513	if (err!=dd_NoError){WerrorS("Error during dd_Matrix2LP in gcone::interiorPoint");}
2514	if (lp==NULL){WerrorS("LP is NULL");}
2515	#ifdef gfan_DEBUG
2516	// dd_WriteLP(stdout,lp);
2517	#endif
2518
2519	lpInt=dd_MakeLPforInteriorFinding(lp);
2520	if (err!=dd_NoError){WerrorS("Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint");}
2521	#ifdef gfan_DEBUG
2522	// dd_WriteLP(stdout,lpInt);
2523	#endif
2524	// dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
2525	if (err!=dd_NoError)
2526	{
2527	WerrorS("Error during dd_FindRelativeInterior in gcone::interiorPoint");
2528	dd_WriteErrorMessages(stdout, err);
2529	}
2530	dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior
2531	// if (err!=dd_NoError){WerrorS("Error during dd_LPSolve");}
2532	lpSol=dd_CopyLPSolution(lpInt);
2533	// if (err!=dd_NoError){WerrorS("Error during dd_CopyLPSolution");}
2534	#ifdef gfan_DEBUG
2535	printf("Interior point: ");
2536	for (int ii=1; ii<(lpSol->d)-1;ii++)
2537	{
2538	dd_WriteNumber(stdout,lpSol->sol[ii]);
2539	}
2540	printf("\n");
2541	#endif
2542	//NOTE The following strongly resembles parts of makeInt.
2543	//Maybe merge sometimes
2544	mpz_t kgV; mpz_init(kgV);
2545	mpz_set_str(kgV,"1",10);
2546	mpz_t den; mpz_init(den);
2547	mpz_t tmp; mpz_init(tmp);
2548	mpq_get_den(tmp,lpSol->sol[1]);
2549	for (int ii=1;ii<(lpSol->d)-1;ii++)
2550	{
2551	mpq_get_den(den,lpSol->sol[ii+1]);
2552	mpz_lcm(kgV,tmp,den);
2553	mpz_set(tmp, kgV);
2554	}
2555	mpq_t qkgV;
2556	mpq_init(qkgV);
2557	mpq_set_z(qkgV,kgV);
2558	for (int ii=1;ii<(lpSol->d)-1;ii++)
2559	{
2560	mpq_t product;
2561	mpq_init(product);
2562	mpq_mul(product,qkgV,lpSol->sol[ii]);
2563	iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
2564	mpq_clear(product);
2565	}
2566	#ifdef gfan_DEBUG
2567	// iv.show();
2568	// cout << endl;
2569	#endif
2570	mpq_clear(qkgV);
2571	mpz_clear(tmp);
2572	mpz_clear(den);
2573	mpz_clear(kgV);
2574
2575	dd_FreeLPSolution(lpSol);
2576	dd_FreeLPData(lpInt);
2577	dd_FreeLPData(lp);
2578	// set_free(ddlinset);
2579	// set_free(ddredrows);
2580
2581	}//void interiorPoint(dd_MatrixPtr const &M)
2582
2583	/** Computes an interior point of a cone by taking two interior points a,b from two different facets
2584	* and then computing b+(a-b)/2
2585	* Of course this only works for flippable facets
2586	* Two cases may occur:
2587	* 1st: There are only two facets who share the only strictly positive ray
2588	* 2nd: There are at least two facets which have a distinct positive ray
2589	* In the former case we use linear algebra to determine an interior point,
2590	* in the latter case we simply add up the two rays
2591	*
2592	* Way too bad! The case may occur that the cone is spanned by three rays, of which only two are strictly
2593	* positive => these lie in a plane and thus their sum is not from relative interior.
2594	* So let's just sum up all rays, find one strictly positive and shift the point along that ray
2595	*
2596	* Used by noRevS
2597	*NOTE no longer used nor maintained. MM Mar 9, 2010
2598	*/
2599	// void gcone::interiorPoint2()
2600	// {//idPrint(this->gcBasis);
2601	// #ifdef gfan_DEBUG
2602	// if(this->ivIntPt!=NULL)
2603	// WarnS("Interior point already exists - ovrewriting!");
2604	// #endif
2605	// facet *f1 = this->facetPtr;
2606	// facet *f2 = NULL;
2607	// intvec *intF1=NULL;
2608	// while(f1!=NULL)
2609	// {
2610	// if(f1->isFlippable)
2611	// {
2612	// facet *f1Ray = f1->codim2Ptr;
2613	// while(f1Ray!=NULL)
2614	// {
2615	// const intvec *check=f1Ray->getRef2FacetNormal();
2616	// if(iv64isStrictlyPositive(check))
2617	// {
2618	// intF1=ivCopy(check);
2619	// break;
2620	// }
2621	// f1Ray=f1Ray->next;
2622	// }
2623	// }
2624	// if(intF1!=NULL)
2625	// break;
2626	// f1=f1->next;
2627	// }
2628	// if(f1!=NULL && f1->next!=NULL)//Choose another facet, different from f1
2629	// f2=f1->next;
2630	// else
2631	// f2=this->facetPtr;
2632	// if(intF1==NULL && hasHomInput==TRUE)
2633	// {
2634	// intF1 = new intvec(this->numVars);
2635	// for(int ii=0;ii<this->numVars;ii++)
2636	// (*intF1)[ii]=1;
2637	// }
2638	// assert(f1); assert(f2);
2639	// intvec *intF2=f2->getInteriorPoint();
2640	// mpq_t *qPosRay = new mpq_t[this->numVars];//The positive ray from above
2641	// mpq_t *qIntPt = new mpq_t[this->numVars];//starting vector a+((b-a)/2)
2642	// mpq_t *qPosIntPt = new mpq_t[this->numVars];//This should be >0 eventually
2643	// for(int ii=0;ii<this->numVars;ii++)
2644	// {
2645	// mpq_init(qPosRay[ii]);
2646	// mpq_init(qIntPt[ii]);
2647	// mpq_init(qPosIntPt[ii]);
2648	// }
2649	// //Compute a+((b-a)/2) && Convert intF1 to mpq
2650	// for(int ii=0;ii<this->numVars;ii++)
2651	// {
2652	// mpq_t a,b;
2653	// mpq_init(a); mpq_init(b);
2654	// mpq_set_si(a,(*intF1)[ii],1);
2655	// mpq_set_si(b,(*intF2)[ii],1);
2656	// mpq_t diff;
2657	// mpq_init(diff);
2658	// mpq_sub(diff,b,a); //diff=b-a
2659	// mpq_t quot;
2660	// mpq_init(quot);
2661	// mpq_div_2exp(quot,diff,1); //quot=diff/2=(b-a)/2
2662	// mpq_clear(diff);
2663	// //Don't be clever and reuse diff here
2664	// mpq_t sum; mpq_init(sum);
2665	// mpq_add(sum,b,quot); //sum=b+quot=a+(b-a)/2
2666	// mpq_set(qIntPt[ii],sum);
2667	// mpq_clear(sum);
2668	// mpq_clear(quot);
2669	// mpq_clear(a); mpq_clear(b);
2670	// //Now for intF1
2671	// mpq_set_si(qPosRay[ii],(*intF1)[ii],1);
2672	// }
2673	// //Now add: qPosIntPt=qPosRay+qIntPt until qPosIntPt >0
2674	// while(TRUE)
2675	// {
2676	// bool success=FALSE;
2677	// int posCtr=0;
2678	// for(int ii=0;ii<this->numVars;ii++)
2679	// {
2680	// mpq_t sum; mpq_init(sum);
2681	// mpq_add(sum,qPosRay[ii],qIntPt[ii]);
2682	// mpq_set(qPosIntPt[ii],sum);
2683	// mpq_clear(sum);
2684	// if(mpq_sgn(qPosIntPt[ii])==1)
2685	// posCtr++;
2686	// }
2687	// if(posCtr==this->numVars)//qPosIntPt > 0
2688	// break;
2689	// else
2690	// {
2691	// mpq_t qTwo; mpq_init(qTwo);
2692	// mpq_set_ui(qTwo,2,1);
2693	// for(int jj=0;jj<this->numVars;jj++)
2694	// {
2695	// mpq_t tmp; mpq_init(tmp);
2696	// mpq_mul(tmp,qPosRay[jj],qTwo);
2697	// mpq_set( qPosRay[jj], tmp);
2698	// mpq_clear(tmp);
2699	// }
2700	// mpq_clear(qTwo);
2701	// }
2702	// }//while
2703	// //Now qPosIntPt ought to be >0, so convert back to int :D
2704	// /Compute lcm of the denominators/
2705	// mpz_t *denom = new mpz_t[this->numVars];
2706	// mpz_t tmp,kgV;
2707	// mpz_init(tmp); mpz_init(kgV);
2708	// for (int ii=0;ii<this->numVars;ii++)
2709	// {
2710	// mpz_t z;
2711	// mpz_init(z);
2712	// mpq_get_den(z,qPosIntPt[ii]);
2713	// mpz_init(denom[ii]);
2714	// mpz_set( denom[ii], z);
2715	// mpz_clear(z);
2716	// }
2717	//
2718	// mpz_set(tmp,denom[0]);
2719	// for (int ii=0;ii<this->numVars;ii++)
2720	// {
2721	// mpz_lcm(kgV,tmp,denom[ii]);
2722	// mpz_set(tmp,kgV);
2723	// }
2724	// mpz_clear(tmp);
2725	// /Multiply the nominators by kgV/
2726	// mpq_t qkgV,res;
2727	// mpq_init(qkgV);
2728	// mpq_canonicalize(qkgV);
2729	// mpq_init(res);
2730	// mpq_canonicalize(res);
2731	//
2732	// mpq_set_num(qkgV,kgV);
2733	// intvec *n=new intvec(this->numVars);
2734	// for (int ii=0;ii<this->numVars;ii++)
2735	// {
2736	// mpq_canonicalize(qPosIntPt[ii]);
2737	// mpq_mul(res,qkgV,qPosIntPt[ii]);
2738	// (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
2739	// }
2740	// this->setIntPoint(n);
2741	// delete n;
2742	// delete [] qPosIntPt;
2743	// delete [] denom;
2744	// delete [] qPosRay;
2745	// delete [] qIntPt;
2746	// mpz_clear(kgV);
2747	// mpq_clear(qkgV); mpq_clear(res);
2748	// }
2749
2750	/** \brief Copy a ring and add a weighted ordering in first place
2751	*
2752	*/
2753	ring gcone::rCopyAndAddWeight(const ring &r, intvec *ivw)
2754	{
2755	ring res=rCopy0(r);
2756	int jj;
2757
2758	omFree(res->order);
2759	res->order =(int )omAlloc0(4sizeof(int));
2760	omFree(res->block0);
2761	res->block0=(int )omAlloc0(4sizeof(int));
2762	omFree(res->block1);
2763	res->block1=(int )omAlloc0(4sizeof(int));
2764	omfree(res->wvhdl);
2765	res->wvhdl =(int *)omAlloc0(4sizeof(int**));
2766
2767	res->order[0]=ringorder_a;
2768	res->block0[0]=1;
2769	res->block1[0]=res->N;;
2770	res->order[1]=ringorder_dp; //basically useless, since that should never be used
2771	res->block0[1]=1;
2772	res->block1[1]=res->N;;
2773	res->order[2]=ringorder_C;
2774
2775	int length=ivw->length();
2776	int A=(int )omAlloc0(length*sizeof(int));
2777	for (jj=0;jj<length;jj++)
2778	{
2779	A[jj]=(*ivw)[jj];
2780	}
2781	res->wvhdl[0]=(int *)A;
2782	res->block1[0]=length;
2783
2784	rComplete(res);
2785	return res;
2786	}//rCopyAndAdd
2787
2788	ring gcone::rCopyAndAddWeight2(const ring &r,const intvec ivw, const intvec fNormal)
2789	{
2790	ring res=rCopy0(r);
2791
2792	omFree(res->order);
2793	res->order =(int )omAlloc0(5sizeof(int));
2794	omFree(res->block0);
2795	res->block0=(int )omAlloc0(5sizeof(int));
2796	omFree(res->block1);
2797	res->block1=(int )omAlloc0(5sizeof(int));
2798	omfree(res->wvhdl);
2799	res->wvhdl =(int *)omAlloc0(5sizeof(int**));
2800
2801	res->order[0]=ringorder_a;
2802	res->block0[0]=1;
2803	res->block1[0]=res->N;
2804	res->order[1]=ringorder_a;
2805	res->block0[1]=1;
2806	res->block1[1]=res->N;
2807
2808	res->order[2]=ringorder_dp;
2809	res->block0[2]=1;
2810	res->block1[2]=res->N;
2811
2812	res->order[3]=ringorder_C;
2813
2814	int length=ivw->length();
2815	int A1=(int )omAlloc0(length*sizeof(int));
2816	int A2=(int )omAlloc0(length*sizeof(int));
2817	for (int jj=0;jj<length;jj++)
2818	{
2819	A1[jj]=(*ivw)[jj];
2820	A2[jj]=-(*fNormal)[jj];
2821	}
2822	res->wvhdl[0]=(int *)A1;
2823	res->block1[0]=length;
2824	res->wvhdl[1]=(int *)A2;
2825	res->block1[1]=length;
2826	rComplete(res);
2827	return res;
2828	}
2829
2830	//NOTE not needed anywhere
2831	// ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
2832	// {
2833	// ring res=rCopy0(currRing);
2834	// rComplete(res);
2835	// rSetWeightVec(res,(int64*)ivw);
2836	// //rChangeCurrRing(rnew);
2837	// return res;
2838	// }
2839
2840	/** \brief Checks whether a given facet is a search facet
2841	* Determines whether a given facet of a cone is the search facet of a neighbouring cone
2842	* This is done in the following way:
2843	* We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
2844	* that is first crossed during the generic walk.
2845	* We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
2846	* If this is the case, then our facet is indeed a search facet and TRUE is retuned.
2847	*/
2848	//removed with r12286
2849
2850	/** \brief Check for equality of two intvecs
2851	*/
2852	static bool ivAreEqual(const intvec &a, const intvec &b)
2853	{
2854	bool res=TRUE;
2855	for(int ii=0;ii<pVariables;ii++)
2856	{
2857	if(a[ii]!=b[ii])
2858	{
2859	res=FALSE;
2860	break;
2861	}
2862	}
2863	return res;
2864	}
2865
2866	/** \brief The reverse search algorithm
2867	*/
2868	//removed with r12286
2869	/** \brief Compute the lineality/homogeneity space
2870	* It is the kernel of the inequality matrix Ax=0
2871	* As a result gcone::dd_LinealitySpace is set
2872	*/
2873	dd_MatrixPtr gcone::computeLinealitySpace()
2874	{
2875	dd_MatrixPtr res;
2876	dd_MatrixPtr ddineq;
2877	ddineq=dd_CopyMatrix(this->ddFacets);
2878	//Add a row of 0s in 0th place
2879	dd_MatrixPtr ddAppendRowOfZeroes=dd_CreateMatrix(1,this->numVars+1);
2880	dd_MatrixPtr ddFoo=dd_AppendMatrix(ddAppendRowOfZeroes,ddineq);
2881	dd_FreeMatrix(ddAppendRowOfZeroes);
2882	dd_FreeMatrix(ddineq);
2883	ddineq=dd_CopyMatrix(ddFoo);
2884	dd_FreeMatrix(ddFoo);
2885	//Cohen starts here
2886	int dimKer=0;//Cohen calls this r
2887	int m=ddineq->rowsize;//Rows
2888	int n=ddineq->colsize;//Cols
2889	int c[m+1];
2890	int d[n+1];
2891	for(int ii=0;ii<m;ii++)
2892	c[ii]=0;
2893	for(int ii=0;ii<n;ii++)
2894	d[ii]=0;
2895
2896	for(int k=1;k<n;k++)
2897	{
2898	//Let's find a j s.t. m[j][k]!=0 && c[j]=0
2899	int condCtr=0;//Check each row for zeroness
2900	for(int j=1;j<m;j++)
2901	{
2902	if(mpq_sgn(ddineq->matrix[j][k])!=0 && c[j]==0)
2903	{
2904	mpq_t quot; mpq_init(quot);
2905	mpq_t one; mpq_init(one); mpq_set_str(one,"-1",10);
2906	mpq_t ratd; mpq_init(ratd);
2907	if((int)mpq_get_d(ddineq->matrix[j][k])!=0)
2908	mpq_div(quot,one,ddineq->matrix[j][k]);
2909	mpq_set(ratd,quot);
2910	mpq_canonicalize(ratd);
2911
2912	mpq_set_str(ddineq->matrix[j][k],"-1",10);
2913	for(int ss=k+1;ss<n;ss++)
2914	{
2915	mpq_t prod; mpq_init(prod);
2916	mpq_mul(prod, ratd, ddineq->matrix[j][ss]);
2917	mpq_set(ddineq->matrix[j][ss],prod);
2918	mpq_canonicalize(ddineq->matrix[j][ss]);
2919	mpq_clear(prod);
2920	}
2921	for(int ii=1;ii<m;ii++)
2922	{
2923	if(ii!=j)
2924	{
2925	mpq_set(ratd,ddineq->matrix[ii][k]);
2926	mpq_set_str(ddineq->matrix[ii][k],"0",10);
2927	for(int ss=k+1;ss<n;ss++)
2928	{
2929	mpq_t prod; mpq_init(prod);
2930	mpq_mul(prod, ratd, ddineq->matrix[j][ss]);
2931	mpq_t sum; mpq_init(sum);
2932	mpq_add(sum, ddineq->matrix[ii][ss], prod);
2933	mpq_set(ddineq->matrix[ii][ss], sum);
2934	mpq_canonicalize(ddineq->matrix[ii][ss]);
2935	mpq_clear(prod);
2936	mpq_clear(sum);
2937	}
2938	}
2939	}
2940	c[j]=k;
2941	d[k]=j;
2942	mpq_clear(quot); mpq_clear(ratd); mpq_clear(one);
2943	}
2944	else
2945	condCtr++;
2946	}
2947	if(condCtr==m-1) //Why -1 ???
2948	{
2949	dimKer++;
2950	d[k]=0;
2951	// break;//goto _4;
2952	}//else{
2953	/Eliminate/
2954	}//for(k
2955	/Output kernel, i.e. set gcone::dd_LinealitySpace here/
2956	// gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
2957	res = dd_CreateMatrix(dimKer,this->numVars+1);
2958	int row=-1;
2959	for(int kk=1;kk<n;kk++)
2960	{
2961	if(d[kk]==0)
2962	{
2963	row++;
2964	for(int ii=1;ii<n;ii++)
2965	{
2966	if(d[ii]>0)
2967	mpq_set(res->matrix[row][ii],ddineq->matrix[d[ii]][kk]);
2968	else if(ii==kk)
2969	mpq_set_str(res->matrix[row][ii],"1",10);
2970	mpq_canonicalize(res->matrix[row][ii]);
2971	}
2972	}
2973	}
2974	dd_FreeMatrix(ddineq);
2975	return(res);
2976	//Better safe than sorry:
2977	//NOTE dd_LinealitySpace contains RATIONAL ENTRIES
2978	//Therefore if you need integer ones: CALL gcone::makeInt(...) method
2979	}
2980
2981
2982	/** \brief The new method of Markwig and Jensen
2983	* Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
2984	* In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
2985	* Keep a list of facets as a linked list containing an intvec and an integer matrix.
2986	* Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
2987	* soon as we compute this facet again. Comparison of facets is done by...
2988	*/
2989	void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
2990	{
2991	facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
2992	// vector<facet> stlRoot;
2993	facet *SearchListAct;
2994	SearchListAct = NULL;
2995	//SearchListAct = SearchListRoot;
2996
2997	gcone *gcAct;
2998	gcAct = &gcRoot;
2999	gcone *gcPtr; //Pointer to end of linked list of cones
3000	gcPtr = &gcRoot;
3001	gcone *gcNext; //Pointer to next cone to be visited
3002	gcNext = &gcRoot;
3003	gcone *gcHead;
3004	gcHead = &gcRoot;
3005
3006	facet *fAct;
3007	fAct = gcAct->facetPtr;
3008
3009	ring rAct;
3010	rAct = currRing;
3011
3012	int UCNcounter=gcAct->getUCN();
3013
3014	#ifdef gfan_DEBUG
3015	printf("NoRevs\n");
3016	printf("Facets are:\n");
3017	gcAct->showFacets();
3018	#endif
3019	/*Compute lineality space here
3020	* Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set*/
3021	if(dd_LinealitySpace==NULL)
3022	dd_LinealitySpace = gcAct->computeLinealitySpace();
3023	/End of lineality space computation/
3024	// gcAct->getCodim2Normals(*gcAct);
3025	if(fAct->codim2Ptr==NULL)
3026	gcAct->getExtremalRays(*gcAct);
3027
3028	//Compute unique representation of Facets and rays, i.e. primitive vectors
3029	// gcAct->normalize();
3030
3031	/* Make a copy of the facet list of first cone
3032	Since the operations getCodim2Normals and normalize affect the facets
3033	we must not memcpy them before these ops! */
3034	/facet codim2Act; codim2Act = NULL;
3035	facet *sl2Root;
3036	facet sl2Act;/
3037	for(int ii=0;ii<this->numFacets;ii++)
3038	{
3039	//only copy flippable facets! NOTE: We do not need flipRing or any such information.
3040	if(fAct->isFlippable==TRUE)
3041	{
3042	/Using shallow copy here/
3043	#ifdef SHALLOW
3044	if( ii==0 \|\| (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
3045	{
3046	if(SearchListRoot!=NULL)
3047	delete(SearchListRoot);
3048	SearchListRoot = fAct->shallowCopy(*fAct);
3049	SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
3050	}
3051	else
3052	{
3053	SearchListAct->next = fAct->shallowCopy(*fAct);
3054	SearchListAct = SearchListAct->next;
3055	}
3056	fAct=fAct->next;
3057	#else
3058	/End of shallow copy/
3059	intvec *fNormal;
3060	fNormal = fAct->getFacetNormal();
3061	if( ii==0 \|\| (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
3062	{
3063	SearchListAct = SearchListRoot; //SearchListRoot is already 'new'ed at beginning of method!
3064	}
3065	else
3066	{
3067	SearchListAct->next = new facet();
3068	SearchListAct = SearchListAct->next;
3069	}
3070	SearchListAct->setFacetNormal(fNormal);
3071	SearchListAct->setUCN(this->getUCN());
3072	SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
3073	SearchListAct->isFlippable=TRUE;
3074	//Copy int point as well
3075	intvec *ivIntPt;// = new intvec(this->numVars);
3076	ivIntPt = fAct->getInteriorPoint();
3077	SearchListAct->setInteriorPoint(ivIntPt);
3078	delete(ivIntPt);
3079	//Copy codim2-facets
3080	facet *codim2Act;
3081	codim2Act = NULL;
3082	facet *sl2Root;
3083	facet *sl2Act;
3084	codim2Act=fAct->codim2Ptr;
3085	SearchListAct->codim2Ptr = new facet(2);
3086	sl2Root = SearchListAct->codim2Ptr;
3087	sl2Act = sl2Root;
3088	//while(codim2Act!=NULL)
3089	for(int jj=0;jj<fAct->numCodim2Facets;jj++)
3090	// for(int jj=0;jj<fAct->numRays-1;jj++)
3091	{
3092	intvec *f2Normal;
3093	f2Normal = codim2Act->getFacetNormal();
3094	if(jj==0)
3095	{
3096	sl2Act = sl2Root;
3097	sl2Act->setFacetNormal(f2Normal);
3098	}
3099	else
3100	{
3101	sl2Act->next = new facet(2);
3102	sl2Act = sl2Act->next;
3103	sl2Act->setFacetNormal(f2Normal);
3104	}
3105	delete f2Normal;
3106	codim2Act = codim2Act->next;
3107	}
3108	fAct = fAct->next;
3109	delete fNormal;
3110	#endif
3111	}//if(fAct->isFlippable==TRUE)
3112	else {fAct = fAct->next;}
3113	}//End of copying facets into SLA
3114
3115	SearchListAct = SearchListRoot; //Set to beginning of list
3116	/Make SearchList doubly linked/
3117	#ifndef NDEBUG
3118	#if SIZEOF_LONG==8
3119	while(SearchListAct!=(facet*)0xfefefefefefefefe && SearchListAct!=NULL)
3120	{
3121	if(SearchListAct->next!=(facet*)0xfefefefefefefefe && SearchListAct->next!=NULL){
3122	#elif SIZEOF_LONG!=8
3123	while(SearchListAct!=(facet*)0xfefefefe)
3124	{
3125	if(SearchListAct->next!=(facet*)0xfefefefe){
3126	#endif
3127	#else
3128	while(SearchListAct!=NULL)
3129	{
3130	if(SearchListAct->next!=NULL){
3131	#endif
3132	// if(SearchListAct->next!=NULL)
3133	// {
3134	SearchListAct->next->prev = SearchListAct;
3135	}
3136	SearchListAct = SearchListAct->next;
3137	}
3138	SearchListAct = SearchListRoot; //Set to beginning of List
3139
3140	// fAct = gcAct->facetPtr;//???
3141	gcAct->writeConeToFile(*gcAct);
3142	/End of initialisation/
3143
3144	fAct = SearchListAct;
3145	/*2nd step
3146	Choose a facet from SearchList, flip it and forget the previous cone
3147	We always choose the first facet from SearchList as facet to be flipped
3148	*/
3149	while( (SearchListAct!=NULL))//&& counter<490)
3150	{//NOTE See to it that the cone is only changed after ALL facets have been flipped!
3151	fAct = SearchListAct;
3152
3153	while(fAct!=NULL)
3154	// while( (fAct->getUCN() == fAct->next->getUCN()) )
3155	{ //Since SLA should only contain flippables there should be no need to check for that
3156	gcAct->flip2(gcAct->gcBasis,fAct);
3157	//NOTE rCopy needed?
3158	ring rTmp=rCopy(fAct->flipRing);
3159	// ring rTmp=fAct->flipRing; //segfaults
3160	rComplete(rTmp);
3161	rChangeCurrRing(rTmp);
3162	gcone gcTmp = new gcone::gcone(gcAct,*fAct);//copy constructor!
3163	/* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
3164	* Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
3165	* NOTE: Can this cause trouble with the destructor?
3166	* Answer: Yes it bloody well does!
3167	* However I'd like to delete this data here, since if we have a cone with many many facets it
3168	* should pay to get rid of the flipGB as soon as possible.
3169	* Destructor must be equipped with necessary checks.
3170	*/
3171	idDelete((ideal *)&fAct->flipGB);
3172	rDelete(fAct->flipRing);
3173
3174	gcTmp->getConeNormals(gcTmp->gcBasis/, FALSE/); //TODO FALSE is default, so should not be needed here
3175	// gcTmp->getCodim2Normals(*gcTmp);
3176	gcTmp->getExtremalRays(*gcTmp);
3177
3178
3179	// //NOTE If flip2 is used we need to get an interior point of gcTmp
3180	// // and replace gcTmp->baseRing with an appropriate ring with only
3181	// // one weight
3182	// gcTmp->interiorPoint2();
3183	gcTmp->replaceDouble_ringorder_a_ByASingleOne();
3184	//gcTmp->ddFacets should not be needed anymore, so
3185	dd_FreeMatrix(gcTmp->ddFacets);
3186	#ifdef gfan_DEBUG
3187	// gcTmp->showFacets(1);
3188	#endif
3189	/add facets to SLA here/
3190	#ifdef SHALLOW
3191	SearchListRoot=gcTmp->enqueue2/NewFacets/(SearchListRoot);
3192	#else
3193	SearchListRoot=gcTmp->enqueueNewFacets(SearchListRoot);
3194	#endif
3195	gcTmp->writeConeToFile(*gcTmp);
3196	if(gfanHeuristic==1)
3197	{
3198	// gcTmp->writeConeToFile(*gcTmp);
3199	idDelete((ideal*)&gcTmp->gcBasis);//Whonder why?
3200	rDelete(gcTmp->baseRing);
3201	}
3202	#ifdef gfan_DEBUG
3203	if(SearchListRoot!=NULL)
3204	showSLA(*SearchListRoot);
3205	#endif
3206	rChangeCurrRing(gcAct->baseRing);
3207	rDelete(rTmp);
3208	//doubly linked for easier removal
3209	gcTmp->prev = gcPtr;
3210	gcPtr->next=gcTmp;
3211	gcPtr=gcPtr->next;
3212	//Cleverly disguised exit condition follows
3213	if(fAct->getUCN() == fAct->next->getUCN())
3214	{
3215	fAct=fAct->next;
3216	}
3217	else
3218	break;
3219	// fAct=fAct->next;
3220	}//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
3221	//Search for cone with smallest UCN
3222	gcNext = gcHead;
3223	#ifndef NDEBUG
3224	#if SIZEOF_LONG==8 //64 bit
3225	while(gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && SearchListRoot!=NULL)
3226	#elif SIZEOF_LONG == 4
3227	while(gcNext!=(gcone * const)0xfbfbfbfb && SearchListRoot!=NULL)
3228	#endif
3229	#endif
3230	#ifdef NDEBUG
3231	while(gcNext!=NULL && SearchListRoot!=NULL)
3232	#endif
3233	{
3234	if( gcNext->getUCN() == SearchListRoot->getUCN() )
3235	{
3236	if(gfanHeuristic==0)
3237	{
3238	gcAct = gcNext;
3239	//Seems better to not to use rCopy here
3240	// rAct=rCopy(gcAct->baseRing);
3241	rAct=gcAct->baseRing;
3242	rComplete(rAct);
3243	rChangeCurrRing(rAct);
3244	break;
3245	}
3246	else if(gfanHeuristic==1)
3247	{
3248	gcone *gcDel;
3249	gcDel = gcAct;
3250	gcAct = gcNext;
3251	//Read st00f from file &
3252	//implant the GB into gcAct
3253	readConeFromFile(gcAct->getUCN(), gcAct);
3254	//Kill the baseRing but ONLY if it is not the ring the computation started in!
3255	// if(gcDel->getUCN()!=1)//WTF!? This causes the Mandelbug in gcone::areEqual(facet, facet)
3256	// rDelete(gcDel->baseRing);
3257	// rAct=rCopy(gcAct->baseRing);
3258	/*The ring change occurs in readConeFromFile, so as to
3259	assure that gcAct->gcBasis belongs to the right ring*/
3260	rAct=gcAct->baseRing;
3261	rComplete(rAct);
3262	rChangeCurrRing(rAct);
3263	break;
3264	}
3265	}
3266	// else if(gcNext->getUCN() < SearchListRoot->getUCN() )
3267	// {
3268	// idDelete( (ideal*)&gcNext->gcBasis );
3269	// rDelete(gcNext->baseRing);
3270	// }
3271	/*else
3272	{
3273	if(gfanHeuristic==1)
3274	{
3275	gcone *gcDel;
3276	gcDel = gcNext;
3277	if(gcDel->getUCN()!=1)
3278	rDelete(gcDel->baseRing);
3279	}
3280	}*/
3281	gcNext = gcNext->next;
3282	}
3283	UCNcounter++;
3284	SearchListAct = SearchListRoot;
3285	}
3286	printf("\nFound %i cones - terminating\n", counter);
3287	}//void noRevS(gcone &gc)
3288
3289
3290	/** \brief Make a set of rational vectors into integers
3291	*
3292	* We compute the lcm of the denominators and multiply with this to get integer values.
3293	* If the gcd of the nominators > 1 we divide by the gcd => primitive vector.
3294	* Expects a new intvec as 3rd parameter
3295	* \param dd_MatrixPtr,intvec
3296	*/
3297	void gcone::makeInt(const dd_MatrixPtr &M, const int line, intvec &n)
3298	{
3299	// mpz_t denom[this->numVars];
3300	mpz_t *denom = new mpz_t[this->numVars];
3301	for(int ii=0;ii<this->numVars;ii++)
3302	{
3303	mpz_init_set_str(denom[ii],"0",10);
3304	}
3305
3306	mpz_t kgV,tmp;
3307	mpz_init(kgV);
3308	mpz_init(tmp);
3309
3310	for (int ii=0;ii<(M->colsize)-1;ii++)
3311	{
3312	mpz_t z;
3313	mpz_init(z);
3314	mpq_get_den(z,M->matrix[line-1][ii+1]);
3315	mpz_set( denom[ii], z);
3316	mpz_clear(z);
3317	}
3318
3319	/Compute lcm of the denominators/
3320	mpz_set(tmp,denom[0]);
3321	for (int ii=0;ii<(M->colsize)-1;ii++)
3322	{
3323	mpz_lcm(kgV,tmp,denom[ii]);
3324	mpz_set(tmp,kgV);
3325	}
3326	mpz_clear(tmp);
3327	/Multiply the nominators by kgV/
3328	mpq_t qkgV,res;
3329	mpq_init(qkgV);
3330	mpq_set_str(qkgV,"1",10);
3331	mpq_canonicalize(qkgV);
3332
3333	mpq_init(res);
3334	mpq_set_str(res,"1",10);
3335	mpq_canonicalize(res);
3336
3337	mpq_set_num(qkgV,kgV);
3338
3339	// mpq_canonicalize(qkgV);
3340	// int ggT=1;
3341	for (int ii=0;ii<(M->colsize)-1;ii++)
3342	{
3343	mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
3344	n[ii]=(int)mpz_get_d(mpq_numref(res));
3345	// ggT=intgcd(ggT,n[ii]);
3346	}
3347	int ggT=n[0];
3348	for(int ii=0;ii<this->numVars;ii++)
3349	ggT=intgcd(ggT,n[ii]);
3350	//Normalization
3351	if(ggT>1)
3352	{
3353	for(int ii=0;ii<this->numVars;ii++)
3354	n[ii] /= ggT;
3355	}
3356	delete [] denom;
3357	mpz_clear(kgV);
3358	mpq_clear(qkgV); mpq_clear(res);
3359
3360	}
3361	/**
3362	* For all codim-2-facets we compute the gcd of the components of the facet normal and
3363	* divide it out. Thus we get a normalized representation of each
3364	* (codim-2)-facet normal, i.e. a primitive vector
3365	* Actually we now also normalize the facet normals.
3366	*/
3367	// void gcone::normalize()
3368	// {
3369	// int *ggT = new int;
3370	// *ggT=1;
3371	// facet *fAct;
3372	// facet *codim2Act;
3373	// fAct = this->facetPtr;
3374	// codim2Act = fAct->codim2Ptr;
3375	// while(fAct!=NULL)
3376	// {
3377	// intvec *fNormal;
3378	// fNormal = fAct->getFacetNormal();
3379	// int *ggT = new int;
3380	// *ggT=1;
3381	// for(int ii=0;ii<this->numVars;ii++)
3382	// {
3383	// ggT=intgcd((ggT),(*fNormal)[ii]);
3384	// }
3385	// if(*ggT>1)//We only need to do this if the ggT is non-trivial
3386	// {
3387	// // intvec *fCopy = fAct->getFacetNormal();
3388	// for(int ii=0;ii<this->numVars;ii++)
3389	// (fNormal)[ii] = ((fNormal)[ii])/(*ggT);
3390	// fAct->setFacetNormal(fNormal);
3391	// }
3392	// delete fNormal;
3393	// delete ggT;
3394	// /And now for the codim2/
3395	// while(codim2Act!=NULL)
3396	// {
3397	// intvec *n;
3398	// n=codim2Act->getFacetNormal();
3399	// int *ggT=new int;
3400	// *ggT=1;
3401	// for(int ii=0;ii<this->numVars;ii++)
3402	// {
3403	// ggT = intgcd((ggT),(*n)[ii]);
3404	// }
3405	// if(*ggT>1)
3406	// {
3407	// for(int ii=0;ii<this->numVars;ii++)
3408	// {
3409	// (n)[ii] = ((n)[ii])/(*ggT);
3410	// }
3411	// codim2Act->setFacetNormal(n);
3412	// }
3413	// codim2Act = codim2Act->next;
3414	// delete n;
3415	// delete ggT;
3416	// }
3417	// fAct = fAct->next;
3418	// }
3419	// }
3420
3421	/** \brief Enqueue new facets into the searchlist
3422	* The searchlist (SLA for short) is implemented as a FIFO
3423	* Checks are done as follows:
3424	* 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
3425	* 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
3426	* facet from SLA and do not add fAct.
3427	* It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
3428	* can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA.
3429	* Takes ptr to search list root
3430	* Returns a pointer to new first element of Searchlist
3431	*/
3432	facet * gcone::enqueueNewFacets(facet *f)
3433	{
3434	#ifdef gfanp
3435	timeval start, end;
3436	gettimeofday(&start, 0);
3437	#endif
3438	facet *slHead;
3439	slHead = f;
3440	facet *slAct; //called with f=SearchListRoot
3441	slAct = f;
3442	facet *slEnd; //Pointer to end of SLA
3443	slEnd = f;
3444	// facet *slEndStatic; //marks the end before a new facet is added
3445	facet *fAct;
3446	fAct = this->facetPtr;
3447	facet *codim2Act;
3448	codim2Act = this->facetPtr->codim2Ptr;
3449	facet *sl2Act;
3450	sl2Act = f->codim2Ptr;
3451	/** Pointer to a facet that will be deleted */
3452	volatile facet *deleteMarker;
3453	deleteMarker = NULL;
3454
3455	/** \brief Flag to mark a facet that might be added
3456	* The following case may occur:
3457	* A facet fAct is found to be parallel but not equal to the current facet slAct from SLA.
3458	* This does however not mean that there does not exist a facet behind the current slAct that is actually equal
3459	* to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
3460	* FALSE and the according slAct is deleted.
3461	* If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
3462	*/
3463	// volatile bool maybe=FALSE;
3464	/**A facet was removed, lengthOfSearchlist-- thus we must not rely on
3465	* if(notParallelCtr==lengthOfSearchList) but rather
3466	* if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
3467	*/
3468	volatile bool removalOccured=FALSE;
3469	// int ctr=0; //encountered equalities in SLA
3470	// int notParallelCtr=0;
3471	// gcone::lengthOfSearchList=1;
3472	while(slEnd->next!=NULL)
3473	{
3474	slEnd=slEnd->next;
3475	// gcone::lengthOfSearchList++;
3476	}
3477	/1st step: compare facetNormals/
3478	while(fAct!=NULL)
3479	{
3480	if(fAct->isFlippable==TRUE)
3481	{
3482	intvec *fNormal=NULL;
3483	fNormal=fAct->getFacetNormal();
3484	slAct = slHead;
3485	//notParallelCtr=0;
3486	/*If slAct==NULL and fAct!=NULL
3487	we just copy all remaining facets into SLA*/
3488	if(slAct==NULL)
3489	{
3490	facet *fCopy;
3491	fCopy = fAct;
3492	while(fCopy!=NULL)
3493	{
3494	if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
3495	{
3496	if(slAct==NULL)
3497	{
3498	slAct = new facet(fCopy/,TRUE*/);//copy constructor
3499	slHead = slAct;
3500	}
3501	else
3502	{
3503	slAct->next = new facet(fCopy/,TRUE*/);
3504	slAct = slAct->next;
3505	}
3506	}
3507	fCopy = fCopy->next;
3508	}
3509	break;//Where does this lead to?
3510	}
3511	/End of dumping into SLA/
3512	while(slAct!=NULL)
3513	{
3514	intvec *slNormal=NULL;
3515	removalOccured=FALSE;
3516	slNormal = slAct->getFacetNormal();
3517	#ifdef gfan_DEBUG
3518	printf("Checking facet (");fNormal->show(1,1);printf(") against (");slNormal->show(1,1);printf(")\n");
3519	#endif
3520	// if( (areEqual(fAct,slAct) && (!areEqual2(fAct,slAct)) ))
3521	// exit(-1);
3522	if(areEqual2(fAct,slAct))
3523	{
3524	deleteMarker = slAct;
3525	if(slAct==slHead)
3526	{
3527	slHead = slAct->next;
3528	if(slHead!=NULL)
3529	slHead->prev = NULL;
3530	}
3531	else if (slAct==slEnd)
3532	{
3533	slEnd=slEnd->prev;
3534	slEnd->next = NULL;
3535	}
3536	else
3537	{
3538	slAct->prev->next = slAct->next;
3539	slAct->next->prev = slAct->prev;
3540	}
3541	removalOccured=TRUE;
3542	gcone::lengthOfSearchList--;
3543	if(deleteMarker!=NULL)
3544	{
3545	// delete deleteMarker;
3546	// deleteMarker=NULL;
3547	}
3548	#ifdef gfan_DEBUG
3549	printf("Removing (");fNormal->show(1,1);printf(") from list\n");
3550	#endif
3551	delete slNormal;
3552	break;//leave the while loop, since we found fAct=slAct thus delete slAct and do not add fAct
3553	}
3554	slAct = slAct->next;
3555	/* NOTE The following lines must not be here but rather called inside the if-block above.
3556	Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now,
3557	(not nice!) but since they are in seperate branches of the if-statement there should not
3558	be a way it gets called twice thus ommiting one facet:
3559	slAct = slAct->next;*/
3560	if(deleteMarker!=NULL)
3561	{
3562	// delete deleteMarker;
3563	// deleteMarker=NULL;
3564	}
3565	delete slNormal;
3566	//if slAct was marked as to be deleted, delete it here!
3567	}//while(slAct!=NULL)
3568	if(removalOccured==FALSE)
3569	{
3570	#ifdef gfan_DEBUG
3571	// cout << "Adding facet (";fNormal->show(1,0);cout << ") to SLA " << endl;
3572	#endif
3573	//Add fAct to SLA
3574	facet *marker;
3575	marker = slEnd;
3576	facet *f2Act;
3577	f2Act = fAct->codim2Ptr;
3578
3579	slEnd->next = new facet();
3580	slEnd = slEnd->next;//Update slEnd
3581	facet *slEndCodim2Root;
3582	facet *slEndCodim2Act;
3583	slEndCodim2Root = slEnd->codim2Ptr;
3584	slEndCodim2Act = slEndCodim2Root;
3585
3586	slEnd->setUCN(this->getUCN());
3587	slEnd->isFlippable = TRUE;
3588	slEnd->setFacetNormal(fNormal);
3589	//NOTE Add interior point here
3590	//This is ugly but needed for flip2
3591	//Better: have slEnd->interiorPoint point to fAct->interiorPoint
3592	//NOTE Only reference -> c.f. copy constructor
3593	//slEnd->setInteriorPoint(fAct->getInteriorPoint());
3594	slEnd->interiorPoint = fAct->interiorPoint;
3595	slEnd->prev = marker;
3596	//Copy codim2-facets
3597	//intvec *f2Normal=new intvec(this->numVars);
3598	while(f2Act!=NULL)
3599	{
3600	intvec *f2Normal;
3601	f2Normal=f2Act->getFacetNormal();
3602	if(slEndCodim2Root==NULL)
3603	{
3604	slEndCodim2Root = new facet(2);
3605	slEnd->codim2Ptr = slEndCodim2Root;
3606	slEndCodim2Root->setFacetNormal(f2Normal);
3607	slEndCodim2Act = slEndCodim2Root;
3608	}
3609	else
3610	{
3611	slEndCodim2Act->next = new facet(2);
3612	slEndCodim2Act = slEndCodim2Act->next;
3613	slEndCodim2Act->setFacetNormal(f2Normal);
3614	}
3615	f2Act = f2Act->next;
3616	delete f2Normal;
3617	}
3618	gcone::lengthOfSearchList++;
3619	}//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) \|\|
3620	fAct = fAct->next;
3621	delete fNormal;
3622	// delete slNormal;
3623	}//if(fAct->isFlippable==TRUE)
3624	else
3625	{
3626	fAct = fAct->next;
3627	}
3628	if(gcone::maxSize<gcone::lengthOfSearchList)
3629	gcone::maxSize= gcone::lengthOfSearchList;
3630	}//while(fAct!=NULL)
3631	#ifdef gfanp
3632	gettimeofday(&end, 0);
3633	time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
3634	#endif
3635	return slHead;
3636	}//addC2N
3637
3638	/** Try using shallow copies*/
3639	facet * gcone::enqueue2(facet *f)
3640	{
3641	assert(f!=NULL);
3642	#ifdef gfanp
3643	timeval start, end;
3644	gettimeofday(&start, 0);
3645	#endif
3646	facet *slHead;
3647	slHead = f;
3648	facet *slAct; //called with f=SearchListRoot
3649	slAct = f;
3650	static facet *slEnd; //Pointer to end of SLA
3651	if(slEnd==NULL)
3652	slEnd = f;
3653
3654	facet *fAct;
3655	fAct = this->facetPtr;//New facets to compare
3656	facet *codim2Act;
3657	codim2Act = this->facetPtr->codim2Ptr;
3658	// facet *sl2Act;
3659	// sl2Act = f->codim2Ptr;
3660	volatile bool removalOccured=FALSE;
3661	while(slEnd->next!=NULL)
3662	{
3663	slEnd=slEnd->next;
3664	}
3665	while(fAct!=NULL)
3666	{
3667	if(fAct->isFlippable)
3668	{
3669	facet *fDeleteMarker=NULL;
3670	slAct = slHead;
3671	if(slAct==NULL)
3672	{
3673	facet *fCopy;
3674	fCopy = fAct;
3675	while(fCopy!=NULL)
3676	{
3677	if(fCopy->isFlippable==TRUE)//We must assure fCopy is also flippable
3678	{
3679	if(slAct==NULL)
3680	{
3681	slAct = slAct->shallowCopy(*fCopy);//shallow copy constructor
3682	slHead = slAct;
3683	}
3684	else
3685	{
3686	slAct->next = slAct->shallowCopy(*fCopy);
3687	slAct = slAct->next;
3688	}
3689	}
3690	fCopy = fCopy->next;
3691	}
3692	break;
3693	}
3694	/Comparison starts here/
3695	while(slAct!=NULL)
3696	{
3697	removalOccured=FALSE;
3698	#ifdef gfan_DEBUG
3699	printf("Checking facet (");fAct->fNormal->show(1,1);printf(") against (");slAct->fNormal->show(1,1);printf(")\n");
3700	#endif
3701	if(areEqual2(fAct,slAct))
3702	{
3703	fDeleteMarker=slAct;
3704	if(slAct==slHead)
3705	{
3706	slHead = slAct->next;
3707	if(slHead!=NULL)
3708	slHead->prev = NULL;
3709	}
3710	else if (slAct==slEnd)
3711	{
3712	slEnd=slEnd->prev;
3713	slEnd->next = NULL;
3714	}
3715	else
3716	{
3717	slAct->prev->next = slAct->next;
3718	slAct->next->prev = slAct->prev;
3719	}
3720	removalOccured=TRUE;
3721	gcone::lengthOfSearchList--;
3722	#ifdef gfan_DEBUG
3723	printf("Removing (");fAct->fNormal->show(1,1);printf(") from list");
3724	#endif
3725	fDeleteMarker->shallowDelete();//Sets everything to NULL
3726	// delete(marker);
3727	break;
3728	}
3729	slAct = slAct->next;
3730	}//while(slAct!=NULL)
3731	if(removalOccured==FALSE)
3732	{
3733	facet *marker=slEnd;
3734	slEnd->next = fAct->shallowCopy(*fAct);
3735	slEnd = slEnd->next;
3736	slEnd->prev=marker;
3737	gcone::lengthOfSearchList++;
3738	}
3739	fAct = fAct->next;
3740	// if(fDeleteMarker!=NULL)
3741	// {
3742	// delete fDeleteMarker;
3743	// fDeleteMarker=NULL;
3744	// }
3745	}
3746	else
3747	fAct = fAct->next;
3748	}//while(fAct!=NULL)
3749
3750	#ifdef gfanp
3751	gettimeofday(&end, 0);
3752	time_enqueue += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
3753	#endif
3754	return slHead;
3755	}
3756
3757	/**
3758	* During flip2 every gc->baseRing gets two ringorder_a
3759	* To avoid having an awful lot of those in the end we endow
3760	* gc->baseRing by a suitable ring with (a,dp,C) and copy all
3761	* necessary stuff over
3762	* But first we will try to just do an inplace exchange and copying only the
3763	* gc->gcBasis
3764	*/
3765	void gcone::replaceDouble_ringorder_a_ByASingleOne()
3766	{
3767	ring srcRing=currRing;
3768	ring replacementRing=rCopy0((ring)this->baseRing);
3769	/We assume we have (a(),a(),dp) here/
3770	omFree(replacementRing->order);
3771	replacementRing->order =(int )omAlloc0(4sizeof(int));
3772	omFree(replacementRing->block0);
3773	replacementRing->block0=(int )omAlloc0(4sizeof(int));
3774	omFree(replacementRing->block1);
3775	replacementRing->block1=(int )omAlloc0(4sizeof(int));
3776	omfree(replacementRing->wvhdl);
3777	replacementRing->wvhdl =(int *)omAlloc0(4sizeof(int**));
3778
3779	replacementRing->order[0]=ringorder_a;
3780	replacementRing->block0[0]=1;
3781	replacementRing->block1[0]=replacementRing->N;
3782
3783	replacementRing->order[1]=ringorder_dp;
3784	replacementRing->block0[1]=1;
3785	replacementRing->block1[1]=replacementRing->N;
3786
3787	replacementRing->order[2]=ringorder_C;
3788
3789	intvec *ivw = this->getIntPoint();
3790	// assert(this->ivIntPt);
3791	int length=ivw->length();
3792	int A=(int )omAlloc0(length*sizeof(int));
3793	for (int jj=0;jj<length;jj++)
3794	A[jj]=(*ivw)[jj];
3795	delete ivw;
3796	replacementRing->wvhdl[0]=(int *)A;
3797	replacementRing->block1[0]=length;
3798	/Finish/
3799	rComplete(replacementRing);
3800	rChangeCurrRing(replacementRing);
3801	ideal temporaryGroebnerBasis=idrCopyR(this->gcBasis,this->baseRing);
3802	rDelete(this->baseRing);
3803	this->baseRing=rCopy(replacementRing);
3804	this->gcBasis=idCopy(temporaryGroebnerBasis);
3805	//FIXME idDelete & rDelete!!! MEMLEAK
3806	/And back to where we came from/
3807	rChangeCurrRing(srcRing);
3808	idDelete( (ideal*)&temporaryGroebnerBasis );
3809	rDelete(replacementRing);
3810	}
3811
3812	/** \brief Compute the gcd of two ints
3813	*/
3814	static int intgcd(const int &a, const int &b)
3815	{
3816	int r, p=a, q=b;
3817	if(p < 0)
3818	p = -p;
3819	if(q < 0)
3820	q = -q;
3821	while(q != 0)
3822	{
3823	r = p % q;
3824	p = q;
3825	q = r;
3826	}
3827	return p;
3828	}
3829
3830	/** \brief Construct a dd_MatrixPtr from a cone's list of facets
3831	* NO LONGER USED
3832	*/
3833	// inline dd_MatrixPtr gcone::facets2Matrix(const gcone &gc)
3834	// {
3835	// facet *fAct;
3836	// fAct = gc.facetPtr;
3837	//
3838	// dd_MatrixPtr M;
3839	// dd_rowrange ddrows;
3840	// dd_colrange ddcols;
3841	// ddcols=(this->numVars)+1;
3842	// ddrows=this->numFacets;
3843	// dd_NumberType numb = dd_Integer;
3844	// M=dd_CreateMatrix(ddrows,ddcols);
3845	//
3846	// int jj=0;
3847	//
3848	// while(fAct!=NULL)
3849	// {
3850	// intvec *comp;
3851	// comp = fAct->getFacetNormal();
3852	// for(int ii=0;ii<this->numVars;ii++)
3853	// {
3854	// dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
3855	// }
3856	// jj++;
3857	// delete comp;
3858	// fAct=fAct->next;
3859	// }
3860	// return M;
3861	// }
3862
3863	/** \brief Write information about a cone into a file on disk
3864	*
3865	* This methods writes the information needed for the "second" method into a file.
3866	* The file's is divided in sections containing information on
3867	* 1) the ring
3868	* 2) the cone's Groebner Basis
3869	* 3) the cone's facets
3870	* Each line contains exactly one date
3871	* Each section starts with its name in CAPITALS
3872	*/
3873	void gcone::writeConeToFile(const gcone &gc, bool usingIntPoints)
3874	{
3875	int UCN=gc.UCN;
3876	stringstream ss;
3877	ss << UCN;
3878	string UCNstr = ss.str();
3879
3880	string prefix="/tmp/Singular/cone";
3881	// string prefix="./"; //crude hack -> run tests in separate directories
3882	string suffix=".sg";
3883	string filename;
3884	filename.append(prefix);
3885	filename.append(UCNstr);
3886	filename.append(suffix);
3887
3888
3889	// int thisPID = getpid();
3890	// ss << thisPID;
3891	// string strPID = ss.str();
3892	// string prefix="./";
3893
3894	ofstream gcOutputFile(filename.c_str());
3895	assert(gcOutputFile);
3896	facet *fAct;
3897	fAct = gc.facetPtr;
3898	facet *f2Act;
3899	f2Act=fAct->codim2Ptr;
3900
3901	char *ringString = rString(gc.baseRing);
3902
3903	if (!gcOutputFile)
3904	{
3905	WerrorS("Error opening file for writing in writeConeToFile\n");
3906	}
3907	else
3908	{
3909	gcOutputFile << "UCN" << endl;
3910	gcOutputFile << this->UCN << endl;
3911	gcOutputFile << "RING" << endl;
3912	gcOutputFile << ringString << endl;
3913	gcOutputFile << "GCBASISLENGTH" << endl;
3914	gcOutputFile << IDELEMS(gc.gcBasis) << endl;
3915	//Write this->gcBasis into file
3916	gcOutputFile << "GCBASIS" << endl;
3917	for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
3918	{
3919	gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
3920	}
3921
3922	gcOutputFile << "FACETS" << endl;
3923
3924	while(fAct!=NULL)
3925	{
3926	const intvec *iv;
3927	iv=fAct->getRef2FacetNormal();//->getFacetNormal();
3928	f2Act=fAct->codim2Ptr;
3929	for (int ii=0;ii<iv->length();ii++)
3930	{
3931	if (ii<iv->length()-1)
3932	{
3933	gcOutputFile << (*iv)[ii] << ",";
3934	}
3935	else
3936	{
3937	gcOutputFile << (*iv)[ii] << " ";
3938	}
3939	}
3940	// delete iv;
3941	while(f2Act!=NULL)
3942	{
3943	const intvec *iv2;
3944	iv2=f2Act->getRef2FacetNormal();//->getFacetNormal();
3945	for(int jj=0;jj<iv2->length();jj++)
3946	{
3947	if (jj<iv2->length()-1)
3948	{
3949	gcOutputFile << (*iv2)[jj] << ",";
3950	}
3951	else
3952	{
3953	gcOutputFile << (*iv2)[jj] << " ";
3954	}
3955	}
3956	// delete iv2;
3957	f2Act = f2Act->next;
3958	}
3959	gcOutputFile << endl;
3960	fAct=fAct->next;
3961	}
3962	gcOutputFile.close();
3963	}
3964
3965	}//writeConeToFile(gcone const &gc)
3966
3967	/** \brief Reads a cone from a file identified by its number
3968	* \|\|depending on whether flip or flip2 is used, switch the flag flipFlag
3969	* \|\|defaults to 0 => flip
3970	* \|\|1 => flip2
3971	*/
3972	void gcone::readConeFromFile(int UCN, gcone *gc)
3973	{
3974	//int UCN=gc.UCN;
3975	stringstream ss;
3976	ss << UCN;
3977	string UCNstr = ss.str();
3978	int gcBasisLength=0;
3979	size_t found; //used for find_first_(not)_of
3980
3981	string prefix="/tmp/Singular/cone";
3982	string suffix=".sg";
3983	string filename;
3984	filename.append(prefix);
3985	filename.append(UCNstr);
3986	filename.append(suffix);
3987
3988	ifstream gcInputFile(filename.c_str(), ifstream::in);
3989
3990	ring saveRing=currRing;
3991	//Comment the following if you uncomment the if(line=="RING") part below
3992	// rChangeCurrRing(gc->baseRing);
3993
3994	while( !gcInputFile.eof() )
3995	{
3996	string line;
3997	getline(gcInputFile,line);
3998	if(line=="RING")
3999	{
4000	getline(gcInputFile,line);
4001	found = line.find("a(");
4002	line.erase(0,found+2);
4003	string strweight;
4004	strweight=line.substr(0,line.find_first_of(")"));
4005
4006	intvec *iv=new intvec(this->numVars);//
4007	for(int ii=0;ii<this->numVars;ii++)
4008	{
4009	string weight;
4010	weight=line.substr(0,line.find_first_of(",)"));
4011	(*iv)[ii]=atol(weight.c_str());//Better to long. Weight bound in Singular:2147483647
4012	line.erase(0,line.find_first_of(",)")+1);
4013	}
4014	found = line.find("a(");
4015
4016	ring newRing;
4017	if(currRing->order[0]!=ringorder_a)
4018	{
4019	newRing=rCopyAndAddWeight(currRing,iv);
4020	}
4021	else
4022	{
4023	newRing=rCopy0(currRing);
4024	int length=this->numVars;
4025	int A=(int )omAlloc0(length*sizeof(int));
4026	for(int jj=0;jj<length;jj++)
4027	{
4028	A[jj]=(*iv)[jj];
4029	}
4030	omFree(newRing->wvhdl[0]);
4031	newRing->wvhdl[0]=(int*)A;
4032	newRing->block1[0]=length;
4033	}
4034	delete iv;
4035	rComplete(newRing);
4036	gc->baseRing=rCopy(newRing);
4037	rDelete(newRing);
4038	rComplete(gc->baseRing);
4039	if(currRing!=gc->baseRing)
4040	rChangeCurrRing(gc->baseRing);
4041	}
4042
4043	if(line=="GCBASISLENGTH")
4044	{
4045	string strGcBasisLength;
4046	getline(gcInputFile, line);
4047	strGcBasisLength = line;
4048	int size=atoi(strGcBasisLength.c_str());
4049	gcBasisLength=size;
4050	gc->gcBasis=idInit(size,1);
4051	}
4052	if(line=="GCBASIS")
4053	{
4054	for(int jj=0;jj<gcBasisLength;jj++)
4055	{
4056	getline(gcInputFile,line);
4057	//magically convert strings into polynomials
4058	//polys.cc:p_Read
4059	//check until first occurance of + or -
4060	//data or c_str
4061	// poly strPoly;//=pInit();//Ought to be inside the while loop, but that will eat your memory
4062	poly resPoly=pInit(); //The poly to be read in
4063	while(!line.empty())
4064	{
4065	poly strPoly;//=pInit();
4066	number nCoeff=nInit(1);
4067	number nCoeffNom=nInit(1);
4068	number nCoeffDenom=nInit(1);
4069	string strMonom, strCoeff, strCoeffNom, strCoeffDenom;
4070	bool hasCoeffInQ = FALSE; //does the polynomial have rational coeff?
4071	bool hasNegCoeff = FALSE; //or a negative one?
4072	found = line.find_first_of("+-"); //get the first monomial
4073	string tmp;
4074	tmp=line[found];
4075	// if(found!=0 && (tmp.compare("-")==0) )
4076	// hasNegCoeff = TRUE; //We have a coeff < 0
4077	if(found==0)
4078	{
4079	if(tmp.compare("-")==0)
4080	hasNegCoeff = TRUE;
4081	line.erase(0,1); //remove leading + or -
4082	found = line.find_first_of("+-"); //adjust found
4083	}
4084	strMonom = line.substr(0,found);
4085	line.erase(0,found);
4086	found = strMonom.find_first_of("/");
4087	if(found!=string::npos) //i.e. "/" exists in strMonom
4088	{
4089	hasCoeffInQ = TRUE;
4090	strCoeffNom=strMonom.substr(0,found);
4091	strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890",found+1));
4092	strMonom.erase(0,found);
4093	strMonom.erase(0,strMonom.find_first_not_of("1234567890/"));
4094	nRead(strCoeffNom.c_str(), &nCoeffNom);
4095	nRead(strCoeffDenom.c_str(), &nCoeffDenom);
4096	}
4097	else
4098	{
4099	found = strMonom.find_first_not_of("1234567890");
4100	strCoeff = strMonom.substr(0,found);
4101	if(!strCoeff.empty())
4102	{
4103	nRead(strCoeff.c_str(),&nCoeff);
4104	}
4105	else
4106	{
4107	// intCoeff = 1;
4108	nCoeff = nInit(1);
4109	}
4110
4111	}
4112	const char* monom = strMonom.c_str();
4113
4114	p_Read(monom,strPoly,currRing); //strPoly:=monom
4115	switch (hasCoeffInQ)
4116	{
4117	case TRUE:
4118	if(hasNegCoeff)
4119	nCoeffNom=nNeg(nCoeffNom);
4120	pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom));
4121	break;
4122	case FALSE:
4123	if(hasNegCoeff)
4124	nCoeff=nNeg(nCoeff);
4125	if(!nIsOne(nCoeff))
4126	{
4127	// if(hasNegCoeff)
4128	// intCoeff *= -1;
4129	// pSetCoeff(strPoly,(number) intCoeff);
4130	pSetCoeff(strPoly, nCoeff );
4131	}
4132	break;
4133
4134	}
4135	//pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1???
4136	if(pIsConstantComp(resPoly))
4137	resPoly=pCopy(strPoly);
4138	else
4139	resPoly=pAdd(resPoly,strPoly);//pAdd = p_Add_q, destroys args
4140	nDelete(&nCoeff);
4141	nDelete(&nCoeffNom);
4142	nDelete(&nCoeffDenom);
4143	// pDelete(&strPoly);
4144	}//while(!line.empty())
4145	gc->gcBasis->m[jj]=pCopy(resPoly);
4146	pDelete(&resPoly); //reset
4147	// pDelete(&strPoly); //NOTE Crashes - already deleted by pAdd
4148	}
4149	// break;
4150	}//if(line=="GCBASIS")
4151	if(line=="FACETS")
4152	{
4153	facet *fAct=gc->facetPtr;
4154	for(int ll=0;ll<this->numFacets;ll++)
4155	{
4156	getline(gcInputFile,line);
4157	found = line.find("\t");
4158	string normalString=line.substr(0,found);
4159	intvec *fN = new intvec(this->numVars);
4160	for(int ii=0;ii<this->numVars;ii++)
4161	{
4162	string component;
4163	found = normalString.find(",");
4164	component=normalString.substr(0,found);
4165	(*fN)[ii]=atol(component.c_str());
4166	normalString.erase(0,found+1);
4167	}
4168	/Only the following line needs to be commented out if you decide not to delete fNormals/
4169	// fAct->setFacetNormal(fN);
4170	delete(fN);
4171	fAct = fAct->next; //Booh, this is ugly
4172	}
4173	break; //NOTE Must always be in the last if-block!
4174	}
4175	}//while(!gcInputFile.eof())
4176	gcInputFile.close();
4177	rChangeCurrRing(saveRing);
4178	}
4179
4180
4181	/** \brief Sort the rays of a facet lexicographically
4182	*/
4183	// void gcone::sortRays(gcone *gc)
4184	// {
4185	// facet *fAct;
4186	// fAct = this->facetPtr->codim2Ptr;
4187	// while(fAct->next!=NULL)
4188	// {
4189	// if(fAct->fNormal->compare(fAct->fNormal->next)==-1
4190	// }
4191	// }
4192
4193	/** \brief Gather the output
4194	* List of lists
4195	* If heuristic==1 readConeFromFile() is called once more on every cone. This may slow down the computation but it also
4196	* allows us to rDelete(gcDel->baseRing) and the such in gcone::noRevS.
4197	\param gc Pointer to gcone, preferably gcRoot ;-)
4198	*\param n the number of cones as determined by gcRoot->getCounter()
4199	*
4200	*/
4201	lists lprepareResult(gcone *gc, const int n)
4202	{
4203	gcone *gcAct;
4204	gcAct = gc;
4205	facet *fAct;
4206	fAct = gc->facetPtr;
4207
4208	lists res=(lists)omAllocBin(slists_bin);
4209	res->Init(n); //initialize to store n cones
4210	for(int ii=0;ii<n;ii++)
4211	{
4212	if(gfanHeuristic==1)// && gcAct->getUCN()>1)
4213	{
4214	gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
4215	// rChangeCurrRing(gcAct->getBaseRing());//NOTE memleak?
4216	}
4217	rChangeCurrRing(gcAct->getRef2BaseRing());
4218	res->m[ii].rtyp=LIST_CMD;
4219	lists l=(lists)omAllocBin(slists_bin);
4220	l->Init(6);
4221	l->m[0].rtyp=INT_CMD;
4222	l->m[0].data=(void*)gcAct->getUCN();
4223	l->m[1].rtyp=IDEAL_CMD;
4224	/*The following is necessary for leaves in the tree of cones
4225	* Since we don't use them in the computation and gcBasis is
4226	* set to (poly)NULL in noRevS we need to get this back here.
4227	*/
4228	// if(gcAct->gcBasis->m[0]==(poly)NULL)
4229	// if(gfanHeuristic==1 && gcAct->getUCN()>1)
4230	// gcAct->readConeFromFile(gcAct->getUCN(),gcAct);
4231	// ring saveRing=currRing;
4232	// ring tmpRing=gcAct->getBaseRing;
4233	// rChangeCurrRing(tmpRing);
4234	// l->m[1].data=(void*)idrCopyR_NoSort(gcAct->gcBasis,gcAct->getBaseRing());
4235	// l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getBaseRing());//NOTE memleak?
4236	l->m[1].data=(void*)idrCopyR(gcAct->gcBasis,gcAct->getRef2BaseRing());
4237	// rChangeCurrRing(saveRing);
4238
4239	l->m[2].rtyp=INTVEC_CMD;
4240	intvec iv=(gcAct->f2M(gcAct,gcAct->facetPtr));//NOTE memleak?
4241	l->m[2].data=(void*)ivCopy(&iv);
4242
4243	l->m[3].rtyp=LIST_CMD;
4244	lists lCodim2List = (lists)omAllocBin(slists_bin);
4245	lCodim2List->Init(gcAct->numFacets);
4246	fAct = gcAct->facetPtr;//fAct->codim2Ptr;
4247	int jj=0;
4248	while(fAct!=NULL && jj<gcAct->numFacets)
4249	{
4250	lCodim2List->m[jj].rtyp=INTVEC_CMD;
4251	intvec ivC2=(gcAct->f2M(gcAct,fAct,2));
4252	lCodim2List->m[jj].data=(void*)ivCopy(&ivC2);
4253	jj++;
4254	fAct = fAct->next;
4255	}
4256	l->m[3].data=(void*)lCodim2List;
4257	l->m[4].rtyp=INTVEC_CMD/RING_CMD/;
4258	l->m[4].data=(void)(gcAct->getIntPoint/BaseRing*/());
4259	l->m[5].rtyp=INT_CMD;
4260	l->m[5].data=(void*)gcAct->getPredUCN();
4261	res->m[ii].data=(void*)l;
4262	gcAct = gcAct->next;
4263	}
4264	return res;
4265	}
4266	/** \brief Write facets of a cone into a matrix
4267	* Takes a pointer to a facet as 2nd arg
4268	* f should always point to gc->facetPtr
4269	* param n is used to determine whether it operates in codim 1 or 2
4270	* We have to cast the int64vecs to intvec due to issues with list structure
4271	*/
4272	inline intvec gcone::f2M(gcone gc, facet f, int n)
4273	{
4274	facet *fAct;
4275	intvec *res;//=new intvec(this->numVars);
4276	// int codim=n;
4277	// int bound;
4278	// if(f==gc->facetPtr)
4279	if(n==1)
4280	{
4281	intvec *m1Res=new intvec(gc->numFacets,gc->numVars,0);
4282	res = ivCopy(m1Res);
4283	fAct = gc->facetPtr;
4284	delete m1Res;
4285	// bound = gc->numFacets*(this->numVars);
4286	}
4287	else
4288	{
4289	fAct = f->codim2Ptr;
4290	intvec *m2Res = new intvec(f->numCodim2Facets,gc->numVars,0);
4291	res = ivCopy(m2Res);
4292	delete m2Res;
4293	// bound = fAct->numCodim2Facets*(this->numVars);
4294
4295	}
4296	int ii=0;
4297	while(fAct!=NULL )//&& ii < bound )
4298	{
4299	const intvec *fNormal;
4300	fNormal = fAct->getRef2FacetNormal();//->getFacetNormal();
4301	for(int jj=0;jj<this->numVars;jj++)
4302	{
4303	(res)[ii]=(int)(fNormal)[jj];//This is ugly and prone to overflow
4304	ii++;
4305	}
4306	fAct = fAct->next;
4307	// delete fNormal;
4308	}
4309	return *res;
4310	}
4311
4312	int gcone::counter=0;
4313	int gfanHeuristic;
4314	int gcone::lengthOfSearchList;
4315	int gcone::maxSize;
4316	dd_MatrixPtr gcone::dd_LinealitySpace;
4317	intvec *gcone::hilbertFunction;
4318	#ifdef gfanp
4319	// int gcone::lengthOfSearchList=0;
4320	float gcone::time_getConeNormals;
4321	float gcone::time_getCodim2Normals;
4322	float gcone::t_getExtremalRays;
4323	float gcone::t_ddPolyh;
4324	float gcone::time_flip;
4325	float gcone::time_flip2;
4326	float gcone::t_areEqual;
4327	float gcone::t_markings;
4328	float gcone::t_dd;
4329	float gcone::t_kStd=0;
4330	float gcone::time_enqueue;
4331	float gcone::time_computeInv;
4332	float gcone::t_ddMC;
4333	float gcone::t_mI;
4334	float gcone::t_iP;
4335	float gcone::t_isParallel;
4336	unsigned gcone::parallelButNotEqual=0;
4337	unsigned gcone::numberOfFacetChecks=0;
4338	#endif
4339	int gcone::numVars;
4340	bool gcone::hasHomInput=FALSE;
4341	intvec *gcone::ivZeroVector;
4342	// ideal gfan(ideal inputIdeal, int h)
4343	lists gfan(ideal inputIdeal, int h)
4344	{
4345	lists lResList; //this is the object we return
4346
4347	if(rHasGlobalOrdering(currRing))
4348	{
4349	// int numvar = pVariables;
4350	gfanHeuristic = h;
4351
4352	enum searchMethod {
4353	reverseSearch,
4354	noRevS
4355	};
4356
4357	searchMethod method;
4358	method = noRevS;
4359
4360	ring inputRing=currRing; // The ring the user entered
4361	// ring rootRing; // The ring associated to the target ordering
4362
4363	dd_set_global_constants();
4364	if(method==noRevS)
4365	{
4366	gcone *gcRoot = new gcone(currRing,inputIdeal);
4367	gcone *gcAct;
4368	gcAct = gcRoot;
4369	gcone::numVars=pVariables;
4370	// gcAct->numVars=pVariables;//NOTE is now static
4371	gcAct->getGB(inputIdeal);
4372	/*Check whether input is homogeneous
4373	if TRUE each facet intersects the positive orthant, so we don't need the
4374	flippability test in getConeNormals & getExtremalRays
4375	*/
4376	if(idHomIdeal(gcAct->gcBasis,NULL))//disabled for tests
4377	{
4378	gcone::hasHomInput=TRUE;
4379	gcone::hilbertFunction=hHstdSeries(inputIdeal,NULL,NULL,NULL,currRing);
4380	}
4381	else
4382	{
4383	gcone::ivZeroVector = new intvec(pVariables);
4384	for(int ii=0;ii<pVariables;ii++)
4385	(*gcone::ivZeroVector)[ii]=0;
4386	}
4387	#ifndef NDEBUG
4388	// cout << "GB of input ideal is:" << endl;
4389	// idShow(gcAct->gcBasis);
4390	#endif
4391	if(isMonomial(gcAct->gcBasis))
4392	{//FIXME
4393	WerrorS("Monomial input - terminating");
4394	dd_free_global_constants();
4395	//This is filthy
4396	goto pointOfNoReturn;
4397	//return lResList;
4398	}
4399	gcAct->getConeNormals(gcAct->gcBasis);
4400	gcone::dd_LinealitySpace = gcAct->computeLinealitySpace();
4401	gcAct->getExtremalRays(*gcAct);
4402	gcAct->noRevS(*gcAct); //Here we go!
4403	//Switch back to the ring the computation was started in
4404	//rChangeCurrRing(inputRing);
4405	//res=gcAct->gcBasis;
4406	//Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
4407	lResList=lprepareResult(gcRoot,gcRoot->getCounter());
4408	/Cleanup/
4409	gcone *gcDel;
4410	gcDel = gcRoot;
4411	gcAct = gcRoot;
4412	while(gcAct!=NULL)
4413	{
4414	gcDel = gcAct;
4415	gcAct = gcAct->next;
4416	// delete gcDel;
4417	}
4418	}//method==noRevS
4419	dd_FreeMatrix(gcone::dd_LinealitySpace);
4420	dd_free_global_constants();
4421	}//rHasGlobalOrdering
4422	else
4423	{
4424	//Simply return an empty list
4425	WerrorS("Ring has non-global ordering.\nThis function requires your current ring to be endowed with a global ordering.\n Now terminating!");
4426	goto pointOfNoReturn;
4427	}
4428	/Return result/
4429	#ifdef gfanp
4430	cout << endl << "t_getConeNormals:" << gcone::time_getConeNormals << endl;
4431	// cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl;
4432	// cout << " t_ddMC:" << gcone::t_ddMC << endl;
4433	// cout << " t_mI:" << gcone::t_mI << endl;
4434	// cout << " t_iP:" << gcone::t_iP << endl;
4435	cout << "t_getExtremalRays:" << gcone::t_getExtremalRays << endl;
4436	cout << " t_ddPolyh:" << gcone::t_ddPolyh << endl;
4437	cout << "t_Flip:" << gcone::time_flip << endl;
4438	cout << " t_markings:" << gcone::t_markings << endl;
4439	cout << " t_dd:" << gcone::t_dd << endl;
4440	cout << " t_kStd:" << gcone::t_kStd << endl;
4441	cout << "t_Flip2:" << gcone::time_flip2 << endl;
4442	cout << " t_dd:" << gcone::t_dd << endl;
4443	cout << " t_kStd:" << gcone::t_kStd << endl;
4444	cout << "t_computeInv:" << gcone::time_computeInv << endl;
4445	cout << "t_enqueue:" << gcone::time_enqueue << endl;
4446	cout << " t_areEqual:" <<gcone::t_areEqual << endl;
4447	cout << "t_isParallel:" <<gcone::t_isParallel << endl;
4448	cout << endl;
4449	cout << "Checked " << gcone::numberOfFacetChecks << " Facets" << endl;
4450	cout << " out of which there were " << gcone::parallelButNotEqual << " parallel but not equal." << endl;
4451	#endif
4452	printf("Maximum lenght of list of facets: %i", gcone::maxSize);
4453	pointOfNoReturn:
4454	return lResList;
4455	}
4456
4457	#endif

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

Download in other formats: