Context Navigation

source: git/kernel/gfan.cc @ 73809b

spielwiese

Last change on this file since 73809b was 73809b, checked in by Martin Monerjan, 14 years ago
#defined SHALLOW facet::shallowCopy, facet::shallowDelete. Needed for gcone::enqueue2 Fixes in gcone::areEqual2(facet,facet), now the standard for comparisons of facet facet::getRef2InteriorPoint getExtremalRays: Fixed dividing out of gcd makeInt: now returns primitve vector renamed gcone::areEqual(intvec&, intvec&) to gcone::ivAreEqual gcone::normalize no longer needed git-svn-id: file:///usr/local/Singular/svn/trunk@12609 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 126.6 KB

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

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

Download in other formats: