Context Navigation

source: git/kernel/gfan.cc @ e58dbb

spielwiese

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

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

Download in other formats: