Context Navigation

source: git/kernel/gfan.cc @ be7ab3f

spielwiese

Last change on this file since be7ab3f was be7ab3f, checked in by Martin Monerjan, 14 years ago
Some asserts Facet copy constructor: no ivCopy of f.interiorPoint, just referencing. Bad style, should go into some kind of shallow copy Skeleton for shallow copy getExtremalRays: For hom input, compute rays, add them up and add a suitable mulitple of (1,_1); for non hom we add the constraints for the pos orthant. Both ways we get strictly pos int points. flip2: Different approach to flipping using an ordering (omega,-N,dp); requires no marking and is appreciably faster than the old way interiorPoint2: Legacy code, no longer used since int points are now computed in getExtremalRays rCopyAndAddWeight2 (a(),a(),dp) Earlier computation of linspace and ext rays for 1st cone interiorPoint goes into the searchlist replaceDouble_ringorder_a_ByASingleOne() changes in readConeFromFile - to be undone git-svn-id: file:///usr/local/Singular/svn/trunk@12589 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 120.1 KB

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

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

Download in other formats: