Context Navigation

source: git/kernel/gfan.cc @ 68eb0c

spielwiese

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

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

Download in other formats: