Context Navigation

source: git/kernel/gfan.cc @ 87d8d5

spielwiese

Last change on this file since 87d8d5 was 87d8d5, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13998 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 123.9 KB

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

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

Download in other formats: