Context Navigation

source: git/kernel/gfan.cc @ 341696

spielwiese

Last change on this file since 341696 was 341696, checked in by Hans Schönemann <hannes@…>, 14 years ago
Adding Id property to all files git-svn-id: file:///usr/local/Singular/svn/trunk@12231 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.4 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009-10-29 16:26:32 $
5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.102 2009-10-29 16:26:32 monerjan Exp $
6	$Id$
7	*/
8
9	#include "mod2.h"
10
11	#ifdef HAVE_GFAN
12
13	#include "kstd1.h"
14	#include "kutil.h"
15	#include "intvec.h"
16	#include "polys.h"
17	#include "ideals.h"
18	#include "kmatrix.h"
19	#include "fast_maps.h" //Mapping of ideals
20	#include "maps.h"
21	#include "ring.h"
22	#include "prCopy.h"
23	#include <iostream>
24	#include <bitset>
25	#include <fstream> //read-write cones to files
26	#include <gmp.h>
27	#include <string>
28	#include <sstream>
29	#include <time.h>
30	//#include <gmpxx.h>
31
32	/DO NOT REMOVE THIS/
33	#ifndef GMPRATIONAL
34	#define GMPRATIONAL
35	#endif
36
37	//Hacks for different working places
38	#define p800
39
40	#ifdef UNI
41	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
42	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
43	#endif
44
45	#ifdef HOME
46	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
47	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
48	#endif
49
50	#ifdef ITWM
51	#include "/u/slg/monerjan/cddlib/include/setoper.h"
52	#include "/u/slg/monerjan/cddlib/include/cdd.h"
53	#include "/u/slg/monerjan/cddlib/include/cddmp.h"
54	#endif
55
56	#ifdef p800
57	#include "../../cddlib/include/setoper.h"
58	#include "../../cddlib/include/cdd.h"
59	#include "../../cddlib/include/cddmp.h"
60	#endif
61
62	#ifndef gfan_DEBUG
63	//#define gfan_DEBUG
64	#ifndef gfan_DEBUGLEVEL
65	#define gfan_DEBUGLEVEL 1
66	#endif
67	#endif
68
69	#include <gfan.h>
70	using namespace std;
71
72	/**
73	*\brief Class facet
74	* Implements the facet structure as a linked list
75	*
76	*/
77
78	/** \brief The default constructor for facets
79	* Do I need a constructor of type facet(intvec)?
80	*/
81	facet::facet()
82	{
83	// Pointer to next facet. */
84	/* Defaults to NULL. This way there is no need to check explicitly */
85	this->fNormal=NULL;
86	this->interiorPoint=NULL;
87	this->UCN=0;
88	this->codim2Ptr=NULL;
89	this->codim=1; //default for (codim-1)-facets
90	this->numCodim2Facets=0;
91	this->flipGB=NULL;
92	this->isIncoming=FALSE;
93	this->next=NULL;
94	this->prev=NULL;
95	this->flipRing=NULL; //the ring on the other side
96	this->isFlippable=FALSE;
97	}
98
99	/** \brief Constructor for facets of codim >= 2
100	*/
101	facet::facet(int const &n)
102	{
103	this->fNormal=NULL;
104	this->interiorPoint=NULL;
105	this->UCN=0;
106	this->codim2Ptr=NULL;
107	if(n>1)
108	{
109	this->codim=n;
110	}//NOTE Handle exception here!
111	this->numCodim2Facets=0;
112	this->flipGB=NULL;
113	this->isIncoming=FALSE;
114	this->next=NULL;
115	this->prev=NULL;
116	this->flipRing=NULL;
117	this->isFlippable=FALSE;
118	}
119
120	/** \brief The copy constructor
121	*/
122	facet::facet(const facet& f)
123	{
124	}
125
126	/** The default destructor */
127	facet::~facet()
128	{
129	delete this->fNormal;
130	delete this->interiorPoint;
131	/* Cleanup the codim2-structure */
132	if(this->codim==2)
133	{
134	facet *codim2Ptr;
135	codim2Ptr = this->codim2Ptr;
136	while(codim2Ptr!=NULL)
137	{
138	delete codim2Ptr->fNormal;
139	codim2Ptr = codim2Ptr->next;
140	}
141	}
142	delete this->codim2Ptr;
143	idDelete((ideal *)&this->flipGB);
144	rDelete(this->flipRing);
145	this->flipRing=NULL;
146	//this=NULL;
147	}
148
149
150	/** \brief Comparison of facets*/
151	bool facet::areEqual(facet &f, facet &g)
152	{
153	bool res = TRUE;
154	intvec *fNormal = new intvec(pVariables);
155	intvec *gNormal = new intvec(pVariables);
156	fNormal = f.getFacetNormal();
157	gNormal = g.getFacetNormal();
158	if((fNormal == gNormal))//\|\|(gcone::isParallel(fNormal,gNormal)))
159	{
160	if(f.numCodim2Facets==g.numCodim2Facets)
161	{
162	facet *f2Act;
163	facet *g2Act;
164	f2Act = f.codim2Ptr;
165	g2Act = g.codim2Ptr;
166	intvec *f2N = new intvec(pVariables);
167	intvec *g2N = new intvec(pVariables);
168	while(f2Act->next!=NULL && g2Act->next!=NULL)
169	{
170	for(int ii=0;ii<f2N->length();ii++)
171	{
172	if(f2Act->getFacetNormal() != g2Act->getFacetNormal())
173	{
174	res=FALSE;
175	}
176	if (res==FALSE)
177	break;
178	}
179	if(res==FALSE)
180	break;
181
182	f2Act = f2Act->next;
183	g2Act = g2Act->next;
184	}//while
185	}//if
186	}//if
187	else
188	{
189	res = FALSE;
190	}
191	return res;
192	}
193
194	/** Stores the facet normal \param intvec*/
195	void facet::setFacetNormal(intvec *iv)
196	{
197	this->fNormal = ivCopy(iv);
198	}
199
200	/** Hopefully returns the facet normal */
201	intvec *facet::getFacetNormal()
202	{
203	return this->fNormal;
204	}
205
206	/** Method to print the facet normal*/
207	void facet::printNormal()
208	{
209	fNormal->show();
210	}
211
212	/** Store the flipped GB*/
213	void facet::setFlipGB(ideal I)
214	{
215	this->flipGB=idCopy(I);
216	}
217
218	/** Return the flipped GB*/
219	ideal facet::getFlipGB()
220	{
221	return this->flipGB;
222	}
223
224	/** Print the flipped GB*/
225	void facet::printFlipGB()
226	{
227	#ifndef NDEBUG
228	idShow(this->flipGB);
229	#endif
230	}
231
232	/** Set the UCN */
233	void facet::setUCN(int n)
234	{
235	this->UCN=n;
236	}
237
238	/** \brief Get the UCN
239	* Returns the UCN iff this != NULL, else -1
240	*/
241	int facet::getUCN()
242	{
243	if(this!=NULL && ( this!=(facet )0xfbfbfbfb \|\| this!=(facet )0xfbfbfbfbfbfbfbfb) )
244	return this->UCN;
245	else
246	return -1;
247	}
248
249	/** Store an interior point of the facet */
250	void facet::setInteriorPoint(intvec *iv)
251	{
252	this->interiorPoint = ivCopy(iv);
253	}
254
255	intvec *facet::getInteriorPoint()
256	{
257	return this->interiorPoint;
258	}
259
260	/** \brief Debugging function
261	* prints the facet normal an all (codim-2)-facets that belong to it
262	*/
263	void facet::fDebugPrint()
264	{
265	facet *codim2Act;
266	codim2Act = this->codim2Ptr;
267	intvec *fNormal;
268	fNormal = this->getFacetNormal();
269	intvec *f2Normal;
270	cout << "=======================" << endl;
271	cout << "Facet normal = (";
272	fNormal->show(1,1);
273	cout << ")"<<endl;
274	cout << "-----------------------" << endl;
275	cout << "Codim2 facets:" << endl;
276	while(codim2Act!=NULL)
277	{
278	f2Normal = codim2Act->getFacetNormal();
279	cout << "(";
280	f2Normal->show(1,0);
281	cout << ")" << endl;
282	codim2Act = codim2Act->next;
283	}
284	cout << "=======================" << endl;
285	}
286
287	//friend class gcone; //Bad style
288
289
290	/**
291	*\brief Implements the cone structure
292	*
293	* A cone is represented by a linked list of facet normals
294	* @see facet
295	*/
296
297
298	/** \brief Default constructor.
299	*
300	* Initialises this->next=NULL and this->facetPtr=NULL
301	*/
302	gcone::gcone()
303	{
304	this->next=NULL;
305	this->prev=NULL;
306	this->facetPtr=NULL; //maybe this->facetPtr = new facet();
307	this->baseRing=currRing;
308	this->counter++;
309	this->UCN=this->counter;
310	this->numFacets=0;
311	this->ivIntPt=NULL;
312	}
313
314	/** \brief Constructor with ring and ideal
315	*
316	* This constructor takes the root ring and the root ideal as parameters and stores
317	* them in the private members gcone::rootRing and gcone::inputIdeal
318	* Since knowledge of the root ring is only needed when using reverse search,
319	* this constructor is not needed when using the "second" method
320	*/
321	gcone::gcone(ring r, ideal I)
322	{
323	this->next=NULL;
324	this->prev=NULL;
325	this->facetPtr=NULL;
326	this->rootRing=r;
327	this->inputIdeal=I;
328	this->baseRing=currRing;
329	this->counter++;
330	this->UCN=this->counter;
331	this->numFacets=0;
332	this->ivIntPt=NULL;
333	}
334
335	/** \brief Copy constructor
336	*
337	* Copies a cone, sets this->gcBasis to the flipped GB
338	* Call this only after a successful call to gcone::flip which sets facet::flipGB
339	*/
340	gcone::gcone(const gcone& gc, const facet &f)
341	{
342	this->next=NULL;
343	// this->prev=(gcone *)&gc; //comment in to get a tree
344	this->prev=NULL;
345	this->numVars=gc.numVars;
346	this->counter++;
347	this->UCN=this->counter;
348	this->facetPtr=NULL;
349	this->gcBasis=idCopy(f.flipGB);
350	this->inputIdeal=idCopy(this->gcBasis);
351	this->baseRing=rCopy(f.flipRing);
352	this->numFacets=0;
353	this->ivIntPt=NULL;
354	this->rootRing=NULL;
355	//rComplete(this->baseRing);
356	//rChangeCurrRing(this->baseRing);
357	}
358
359	/** \brief Default destructor
360	*/
361	gcone::~gcone()
362	{
363	if(this->gcBasis!=NULL)
364	idDelete((ideal *)&this->gcBasis);
365	if(this->inputIdeal!=NULL)
366	idDelete((ideal *)&this->inputIdeal);
367	if (this->rootRing!=NULL && this->rootRing!=(ip_sring *)0xfefefefefefefefe)
368	rDelete(this->rootRing);
369	//rDelete(this->baseRing);
370	facet *fAct;
371	facet *fDel;
372	/Delete the facet structure/
373	fAct=this->facetPtr;
374	fDel=fAct;
375	while(fAct!=NULL)
376	{
377	fDel=fAct;
378	fAct=fAct->next;
379	delete fDel;
380	}
381	}
382
383	/** \brief Set the interior point of a cone */
384	void gcone::setIntPoint(intvec *iv)
385	{
386	this->ivIntPt=ivCopy(iv);
387	}
388
389	/** \brief Return the interior point */
390	intvec *gcone::getIntPoint()
391	{
392	return this->ivIntPt;
393	}
394
395	/** \brief Print the interior point */
396	void gcone::showIntPoint()
397	{
398	ivIntPt->show();
399	}
400
401	/** \brief Print facets
402	* This is mainly for debugging purposes. Usually called from within gdb
403	*/
404	void gcone::showFacets(short codim)
405	{
406	facet *f;
407	switch(codim)
408	{
409	case 1:
410	f = this->facetPtr;
411	break;
412	case 2:
413	f = this->facetPtr->codim2Ptr;
414	break;
415	}
416
417	intvec *iv;
418	while(f!=NULL)
419	{
420	iv = f->getFacetNormal();
421	cout << "(";
422	iv->show(1,0);
423	if(f->isFlippable==FALSE)
424	cout << ")* ";
425	else
426	cout << ") ";
427	f=f->next;
428	}
429	cout << endl;
430	}
431
432	/** For debugging purposes only */
433	void gcone::showSLA(facet &f)
434	{
435	facet *fAct;
436	fAct = &f;
437	facet *codim2Act;
438	codim2Act = fAct->codim2Ptr;
439	intvec *fNormal;
440	intvec *f2Normal;
441	cout << endl;
442	while(fAct!=NULL)
443	{
444	fNormal=fAct->getFacetNormal();
445	cout << "(";
446	fNormal->show(1,0);
447	if(fAct->isFlippable==TRUE)
448	cout << ") ";
449	else
450	cout << ")* ";
451	codim2Act = fAct->codim2Ptr;
452	cout << " Codim2: ";
453	while(codim2Act!=NULL)
454	{
455	f2Normal = codim2Act->getFacetNormal();
456	cout << "(";
457	f2Normal->show(1,0);
458	cout << ") ";
459	codim2Act = codim2Act->next;
460	}
461	cout << "UCN = " << fAct->getUCN() << endl;
462	fAct = fAct->next;
463	}
464	}
465
466	void gcone::idDebugPrint(ideal const &I)
467	{
468	int numElts=IDELEMS(I);
469	cout << "Ideal with " << numElts << " generators" << endl;
470	cout << "Leading terms: ";
471	for (int ii=0;ii<numElts;ii++)
472	{
473	pWrite0(pHead(I->m[ii]));
474	cout << ",";
475	}
476	cout << endl;
477	}
478
479	bool gcone::isMonomial(ideal const &I)
480	{
481	bool res = TRUE;
482	for(int ii=0;ii<IDELEMS(I);ii++)
483	{
484	if(pLength((poly)I->m[ii])>1)
485	{
486	res = FALSE;
487	break;
488	}
489	}
490	return res;
491	}
492
493	/** \brief Set gcone::numFacets */
494	void gcone::setNumFacets()
495	{
496	}
497
498	/** \brief Get gcone::numFacets */
499	int gcone::getNumFacets()
500	{
501	return this->numFacets;
502	}
503
504	int gcone::getUCN()
505	{
506	if( this!=NULL && ( this!=(gcone * const)0xfbfbfbfbfbfbfbfb && this!=(gcone * const)0xfbfbfbfb ) )
507	return this->UCN;
508	else
509	return -1;
510	}
511	/** \brief Compute the normals of the cone
512	*
513	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
514	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
515	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
516	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
517	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
518	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
519	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
520	*
521	* Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
522	* an interior point of the cone.
523	*/
524	void gcone::getConeNormals(ideal const &I, bool compIntPoint)
525	{
526	#ifdef gfan_DEBUG
527	std::cout << "* Computing Inequalities... *" << std::endl;
528	#endif
529	//All variables go here - except ineq matrix and *v, see below
530	int lengthGB=IDELEMS(I); // # of polys in the groebner basis
531	int pCompCount; // # of terms in a poly
532	poly aktpoly;
533	int numvar = pVariables; // # of variables in a polynomial (or ring?)
534	int leadexp[numvar]; // dirty hack of exp.vects
535	int aktexp[numvar];
536	int cols,rows; // will contain the dimensions of the ineq matrix - deprecated by
537	dd_rowrange ddrows;
538	dd_colrange ddcols;
539	dd_rowset ddredrows; // # of redundant rows in ddineq
540	dd_rowset ddlinset; // the opposite
541	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
542	dd_NumberType ddnumb=dd_Integer; //Number type
543	dd_ErrorType dderr=dd_NoError; //
544	// End of var declaration
545	#ifdef gfan_DEBUG
546	// cout << "The Groebner basis has " << lengthGB << " elements" << endl;
547	// cout << "The current ring has " << numvar << " variables" << endl;
548	#endif
549	cols = numvar;
550
551	//Compute the # inequalities i.e. rows of the matrix
552	rows=0; //Initialization
553	for (int ii=0;ii<IDELEMS(I);ii++)
554	{
555	aktpoly=(poly)I->m[ii];
556	rows=rows+pLength(aktpoly)-1;
557	}
558	#ifdef gfan_DEBUG
559	// cout << "rows=" << rows << endl;
560	// cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
561	#endif
562	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
563	dd_set_global_constants();
564	ddrows=rows;
565	ddcols=cols;
566	dd_MatrixPtr ddineq; //Matrix to store the inequalities
567	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
568
569	// We loop through each g\in GB and compute the resulting inequalities
570	for (int i=0; i<IDELEMS(I); i++)
571	{
572	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
573	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
574	#ifdef gfan_DEBUG
575	// cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
576	#endif
577
578	int v=(int )omAlloc((numvar+1)*sizeof(int));
579	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
580
581	//Store leadexp for aktpoly
582	for (int kk=0;kk<numvar;kk++)
583	{
584	leadexp[kk]=v[kk+1];
585	//Since we need to know the difference of leadexp with the other expvects we do nothing here
586	//but compute the diff below
587	}
588
589
590	while (pNext(aktpoly)!=NULL) //move to next term until NULL
591	{
592	aktpoly=pNext(aktpoly);
593	pSetm(aktpoly); //doesn't seem to help anything
594	pGetExpV(aktpoly,v);
595	for (int kk=0;kk<numvar;kk++)
596	{
597	aktexp[kk]=v[kk+1];
598	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
599	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
600	}
601	aktmatrixrow=aktmatrixrow+1;
602	}//while
603
604	} //for
605
606	//Maybe add another row to contain the constraints of the standard simplex?
607
608	#ifdef gfan_DEBUG
609	// cout << "The inequality matrix is" << endl;
610	// dd_WriteMatrix(stdout, ddineq);
611	#endif
612
613	// The inequalities are now stored in ddineq
614	// Next we check for superflous rows
615	ddredrows = dd_RedundantRows(ddineq, &dderr);
616	if (dderr!=dd_NoError) // did an error occur?
617	{
618	dd_WriteErrorMessages(stderr,dderr); //if so tell us
619	}
620	#ifdef gfan_DEBUG
621	else
622	{
623	// cout << "Redundant rows: ";
624	// set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
625	}//if dd_Error
626	#endif
627	//Remove reduntant rows here!
628	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
629	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
630	ddcols = ddineq->colsize;
631
632	//ddCreateMatrix(ddrows,ddcols+1);
633	ddFacets = dd_CopyMatrix(ddineq);
634	#ifdef gfan_DEBUG
635	// cout << "Having removed redundancies, the normals now read:" << endl;
636	// dd_WriteMatrix(stdout,ddineq);
637	// cout << "Rows = " << ddrows << endl;
638	// cout << "Cols = " << ddcols << endl;
639	#endif
640
641	/Write the normals into class facet/
642	#ifdef gfan_DEBUG
643	// cout << "Creating list of normals" << endl;
644	#endif
645	/The pointer fRoot should be the return value of this function*/
646	//facet *fRoot = new facet(); //instantiate new facet
647	//fRoot->setUCN(this->getUCN());
648	//this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
649	facet *fAct; //pointer to active facet
650	//fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress!
651	//std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
652	int numNonFlip=0;
653	for (int kk = 0; kk<ddrows; kk++)
654	{
655	intvec *load = new intvec(this->numVars); //intvec to store a single facet normal that will then be stored via setFacetNormal
656	for (int jj = 1; jj <ddcols; jj++)
657	{
658	double foo;
659	foo = mpq_get_d(ddineq->matrix[kk][jj]);
660	#ifdef gfan_DEBUG
661	// std::cout << "fAct is " << foo << " at " << fAct << std::endl;
662	#endif
663	(*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
664	}//for (int jj = 1; jj <ddcols; jj++)
665
666	/Quick'n'dirty hack for flippability/
667	bool isFlip=FALSE;
668
669	for (int jj = 0; jj<load->length(); jj++)
670	{
671	intvec *ivCanonical = new intvec(load->length());
672	(*ivCanonical)[jj]=1;
673	// cout << "dotProd=" << dotProduct(load,ivCanonical) << endl;
674	if (dotProduct(load,ivCanonical)<0)
675	//if (ivMult(load,ivCanonical)<0)
676	{
677	isFlip=TRUE;
678	break; //URGHS
679	}
680	}
681	if (isFlip==FALSE)
682	{
683	this->numFacets++;
684	numNonFlip++;
685	if(this->numFacets==1)
686	{
687	facet *fRoot = new facet();
688	this->facetPtr = fRoot;
689	fAct = fRoot;
690
691	}
692	else
693	{
694	fAct->next = new facet();
695	fAct = fAct->next;
696	}
697	fAct->isFlippable=FALSE;
698	fAct->setFacetNormal(load);
699	fAct->setUCN(this->getUCN());
700	#ifdef gfan_DEBUG
701	cout << "Marking facet (";
702	load->show(1,0);
703	cout << ") as non flippable" << endl;
704	#endif
705	}
706	else
707	{ /Now load should be full and we can call setFacetNormal/
708	this->numFacets++;
709	if(this->numFacets==1)
710	{
711	facet *fRoot = new facet();
712	this->facetPtr = fRoot;
713	fAct = fRoot;
714	}
715	else
716	{
717	fAct->next = new facet();
718	fAct = fAct->next;
719	}
720	fAct->isFlippable=TRUE;
721	fAct->setFacetNormal(load);
722	fAct->setUCN(this->getUCN());
723	}//if (isFlippable==FALSE)
724	//delete load;
725	}//for (int kk = 0; kk<ddrows; kk++)
726
727	//In cases like I=<x-1,y-1> there are only non-flippable facets...
728	if(numNonFlip==this->numFacets)
729	{
730	cerr << "Only non-flippable facets. Terminating..." << endl;
731	}
732
733	/*
734	Now we should have a linked list containing the facet normals of those facets that are
735	-irredundant
736	-flipable
737	Adressing is done via *facetPtr
738	*/
739
740	if (compIntPoint==TRUE)
741	{
742	intvec *iv = new intvec(this->numVars);
743	dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
744	int jj=1;
745	for (int ii=0;ii<=this->numVars;ii++)
746	{
747	dd_set_si(posRestr->matrix[ii][jj],1);
748	jj++;
749	}
750	dd_MatrixAppendTo(&ddineq,posRestr);
751	interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets
752	this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt
753	//delete iv;
754	}
755
756	//Compute the number of facets
757	// wrong because ddineq->rowsize contains only irredundant facets. There can still be
758	// non-flippable ones!
759	//this->numFacets=ddineq->rowsize;
760
761	//Clean up but don't delete the return value! (Whatever it will turn out to be)
762	dd_FreeMatrix(ddineq);
763	set_free(ddredrows);
764	//free(ddnewpos);
765	//set_free(ddlinset);
766	//NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
767	//THIS SUCKS
768	//dd_free_global_constants();
769
770	}//method getConeNormals(ideal I)
771
772	/** \brief Compute the (codim-2)-facets of a given cone
773	* This method is used during noRevS
774	* Additionally we check whether the codim2-facet normal is strictly positive. Otherwise
775	* the facet is marked as non-flippable.
776	*/
777	void gcone::getCodim2Normals(gcone const &gc)
778	{
779	//this->facetPtr->codim2Ptr = new facet(2); //instantiate a (codim-2)-facet
780	facet *fAct;
781	fAct = this->facetPtr;
782	facet *codim2Act;
783	//codim2Act = this->facetPtr->codim2Ptr;
784
785	dd_set_global_constants();
786
787	dd_MatrixPtr ddineq,P,ddakt;
788	dd_rowset impl_linset, redset;
789	dd_ErrorType err;
790	dd_rowindex newpos;
791
792	//ddineq = facets2Matrix(gc); //get a matrix representation of the cone
793	ddineq = dd_CopyMatrix(gc.ddFacets);
794
795	/Now set appropriate linearity/
796	dd_PolyhedraPtr ddpolyh;
797	for (int ii=0; ii<this->numFacets; ii++)
798	{
799	ddakt = dd_CopyMatrix(ddineq);
800	ddakt->representation=dd_Inequality;
801	set_addelem(ddakt->linset,ii+1);
802	dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
803
804	#ifdef gfan_DEBUG
805	// cout << "Codim2 matrix"<<endl;
806	// dd_WriteMatrix(stdout,ddakt);
807	#endif
808	ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
809	P=dd_CopyGenerators(ddpolyh);
810	#ifdef gfan_DEBUG
811	// cout << "Codim2 facet:" << endl;
812	// dd_WriteMatrix(stdout,P);
813	// cout << endl;
814	#endif
815
816	/* We loop through each row of P
817	* normalize it by making all entries integer ones
818	* and add the resulting vector to the int matrix facet::codim2Facets
819	*/
820	for (int jj=1;jj<=P->rowsize;jj++)
821	{
822	fAct->numCodim2Facets++;
823	if(fAct->numCodim2Facets==1)
824	{
825	fAct->codim2Ptr = new facet(2);
826	codim2Act = fAct->codim2Ptr;
827	}
828	else
829	{
830	codim2Act->next = new facet(2);
831	codim2Act = codim2Act->next;
832	}
833	intvec *n = new intvec(this->numVars);
834	makeInt(P,jj,*n);
835	codim2Act->setFacetNormal(n);
836	delete n;
837	}
838	/*We check whether the facet spanned by the codim-2 facets
839	* intersects with the positive orthant. Otherwise we define this
840	* facet to be non-flippable
841	*/
842	intvec *iv_intPoint = new intvec(this->numVars);
843	dd_MatrixPtr shiftMatrix;
844	dd_MatrixPtr intPointMatrix;
845	shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
846	for(int kk=0;kk<this->numVars;kk++)
847	{
848	dd_set_si(shiftMatrix->matrix[kk][0],1);
849	dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
850	}
851	intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
852	//dd_WriteMatrix(stdout,intPointMatrix);
853	interiorPoint(intPointMatrix,*iv_intPoint);
854	for(int ll=0;ll<this->numVars;ll++)
855	{
856	if( (*iv_intPoint)[ll] < 0 )
857	{
858	fAct->isFlippable=FALSE;
859	break;
860	}
861	}
862	/End of check/
863	fAct = fAct->next;
864	dd_FreeMatrix(ddakt);
865	// dd_FreeMatrix(ddineq);
866	dd_FreePolyhedra(ddpolyh);
867	delete iv_intPoint;
868	}//while
869	}
870
871	/** \brief Compute the Groebner Basis on the other side of a shared facet
872	*
873	* Implements algorithm 4.3.2 from Anders' thesis.
874	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
875	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
876	* is parallel to \f$ leadexp(g) \f$
877	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
878	* Other possibilities include computing the rank of the matrix consisting of the vectors in question and
879	* computing an interior point of the facet and taking all terms having the same weight with respect
880	* to this interior point.
881	*\param ideal, facet
882	* Input: a marked,reduced Groebner basis and a facet
883	*/
884	void gcone::flip(ideal gb, facet *f) //Compute "the other side"
885	{
886	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
887	fNormal = f->getFacetNormal(); //read this->fNormal;
888	//#ifdef gfan_DEBUG
889	//std::cout << "===" << std::endl;
890	std::cout << "running gcone::flip" << std::endl;
891	std::cout << "flipping UCN " << this->getUCN() << endl;
892	// for(int ii=0;ii<IDELEMS(gb);ii++) //not very handy with large examples
893	// {
894	// pWrite((poly)gb->m[ii]);
895	// }
896	cout << "over facet (";
897	fNormal->show(1,0);
898	cout << ") with UCN " << f->getUCN();
899	std::cout << std::endl;
900	//#endif
901	/1st step: Compute the initial ideal/
902	poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form
903	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
904	poly aktpoly;
905	intvec *check = new intvec(this->numVars); //array to store the difference of LE and v
906
907	for (int ii=0;ii<IDELEMS(gb);ii++)
908	{
909	aktpoly = (poly)gb->m[ii];
910	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
911	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
912	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
913	initialFormElement[ii]=pHead(aktpoly);
914
915	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
916	{
917	aktpoly=pNext(aktpoly); //next term
918	pSetm(aktpoly);
919	pGetExpV(aktpoly,v);
920	/* Convert (int)v into (intvec)check */
921	for (int jj=0;jj<this->numVars;jj++)
922	{
923	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
924	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
925	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
926	}
927	#ifdef gfan_DEBUG
928	// cout << "check=";
929	// check->show();
930	// cout << endl;
931	#endif
932	if (isParallel(check,fNormal)) //pass *check when
933	{
934	// cout << "Parallel vector found, adding to initialFormElement" << endl;
935	initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
936	}
937	}//while
938	#ifdef gfan_DEBUG
939	// cout << "Initial Form=";
940	// pWrite(initialFormElement[ii]);
941	// cout << "---" << endl;
942	#endif
943	/Now initialFormElement must be added to (ideal)initialForm /
944	initialForm->m[ii]=initialFormElement[ii];
945	}//for
946	#ifdef gfan_DEBUG
947	/* cout << "Initial ideal is: " << endl;
948	idShow(initialForm);
949	//f->printFlipGB();*/
950	// cout << "===" << endl;
951	#endif
952	//delete check;
953
954	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
955	/*Substep 2.1
956	compute $G_{-\alpha}(in_v(I))
957	see journal p. 66
958	NOTE Check for different rings. Prolly it will not always be necessary to add a weight, if the
959	srcRing already has a weighted ordering
960	*/
961	ring srcRing=currRing;
962	ring tmpRing;
963
964	if( (srcRing->order[0]!=ringorder_a))
965	{
966	tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
967	}
968	else
969	{
970	tmpRing=rCopy0(srcRing);
971	int length=fNormal->length();
972	int A=(int )omAlloc0(length*sizeof(int));
973	for(int jj=0;jj<length;jj++)
974	{
975	A[jj]=-(*fNormal)[jj];
976	}
977	omFree(tmpRing->wvhdl[0]);
978	tmpRing->wvhdl[0]=(int*)A;
979	tmpRing->block1[0]=length;
980	rComplete(tmpRing);
981	}
982	rChangeCurrRing(tmpRing);
983
984	//rWrite(currRing); cout << endl;
985
986	ideal ina;
987	ina=idrCopyR(initialForm,srcRing);
988	#ifdef gfan_DEBUG
989	// cout << "ina=";
990	// idShow(ina); cout << endl;
991	#endif
992	ideal H;
993	//H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
994	H=kStd(ina,NULL,testHomog,NULL);
995	idSkipZeroes(H);
996	#ifdef gfan_DEBUG
997	// cout << "H="; idShow(H); cout << endl;
998	#endif
999	/*Substep 2.2
1000	do the lifting and mark according to H
1001	*/
1002	rChangeCurrRing(srcRing);
1003	ideal srcRing_H;
1004	ideal srcRing_HH;
1005	srcRing_H=idrCopyR(H,tmpRing);
1006	#ifdef gfan_DEBUG
1007	// cout << "srcRing_H = ";
1008	// idShow(srcRing_H); cout << endl;
1009	#endif
1010	srcRing_HH=ffG(srcRing_H,this->gcBasis);
1011	#ifdef gfan_DEBUG
1012	// cout << "srcRing_HH = ";
1013	// idShow(srcRing_HH); cout << endl;
1014	#endif
1015	/*Substep 2.2.1
1016	* Mark according to G_-\alpha
1017	* Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
1018	* we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
1019	* represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
1020	* Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
1021	* compute the difference accordingly
1022	*/
1023	dd_set_global_constants();
1024	bool markingsAreCorrect=FALSE;
1025	dd_MatrixPtr intPointMatrix;
1026	int iPMatrixRows=0;
1027	dd_rowrange aktrow=0;
1028	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1029	{
1030	poly aktpoly=(poly)srcRing_HH->m[ii];
1031	iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
1032	}
1033	/* additionally one row for the standard-simplex and another for a row that becomes 0 during
1034	* construction of the differences
1035	*/
1036	intPointMatrix = dd_CreateMatrix(iPMatrixRows+10,this->numVars+1); //iPMatrixRows+10;
1037	intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
1038
1039	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1040	{
1041	markingsAreCorrect=FALSE; //crucial to initialise here
1042	poly aktpoly=srcRing_HH->m[ii];
1043	/Comparison of leading monomials is done via exponent vectors/
1044	for (int jj=0;jj<IDELEMS(H);jj++)
1045	{
1046	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
1047	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
1048	pGetExpV(aktpoly,src_ExpV);
1049	rChangeCurrRing(tmpRing); //this ring change is crucial!
1050	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
1051	rChangeCurrRing(srcRing);
1052	bool expVAreEqual=TRUE;
1053	for (int kk=1;kk<=this->numVars;kk++)
1054	{
1055	#ifdef gfan_DEBUG
1056	cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
1057	#endif
1058	if (src_ExpV[kk]!=dst_ExpV[kk])
1059	{
1060	expVAreEqual=FALSE;
1061	}
1062	}
1063	//if (src_ExpV == dst_ExpV)
1064	if (expVAreEqual==TRUE)
1065	{
1066	markingsAreCorrect=TRUE; //everything is fine
1067	#ifdef gfan_DEBUG
1068	cout << "correct markings" << endl;
1069	#endif
1070	}//if (pHead(aktpoly)==pHead(H->m[jj])
1071	delete src_ExpV;
1072	delete dst_ExpV;
1073	}//for (int jj=0;jj<IDELEMS(H);jj++)
1074
1075	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1076	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
1077	if (markingsAreCorrect==TRUE)
1078	{
1079	pGetExpV(aktpoly,leadExpV);
1080	}
1081	else
1082	{
1083	rChangeCurrRing(tmpRing);
1084	pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
1085	rChangeCurrRing(srcRing);
1086	}
1087	/compute differences of the expvects/
1088	while (pNext(aktpoly)!=NULL)
1089	{
1090	/*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
1091	is not omitted when computing the differences*/
1092	if(markingsAreCorrect==TRUE)
1093	{
1094	aktpoly=pNext(aktpoly);
1095	pGetExpV(aktpoly,v);
1096	}
1097	else
1098	{
1099	pGetExpV(pHead(aktpoly),v);
1100	markingsAreCorrect=TRUE;
1101	}
1102
1103	for (int jj=0;jj<this->numVars;jj++)
1104	{
1105	/Store into ddMatrix/
1106	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
1107	}
1108	aktrow +=1;
1109	}
1110	delete v;
1111	delete leadExpV;
1112	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
1113	/Now we add the constraint for the standard simplex/
1114	/NOTE:Might actually work without the standard simplex/
1115	dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
1116	for (int jj=1;jj<=this->numVars;jj++)
1117	{
1118	dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
1119	}
1120	//Let's make sure we compute interior points from the positive orthant
1121	dd_MatrixPtr posRestr=dd_CreateMatrix(this->numVars,this->numVars+1);
1122	int jj=1;
1123	for (int ii=0;ii<this->numVars;ii++)
1124	{
1125	dd_set_si(posRestr->matrix[ii][jj],1);
1126	jj++;
1127	}
1128	dd_MatrixAppendTo(&intPointMatrix,posRestr);
1129	dd_FreeMatrix(posRestr);
1130	#ifdef gfan_DEBUG
1131	// dd_WriteMatrix(stdout,intPointMatrix);
1132	#endif
1133	intvec *iv_weight = new intvec(this->numVars);
1134	interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point
1135	dd_FreeMatrix(intPointMatrix);
1136	dd_free_global_constants();
1137
1138	/*Step 3
1139	turn the minimal basis into a reduced one
1140	*/
1141	//ring dstRing=rCopyAndAddWeight(tmpRing,iv_weight);
1142
1143	// NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
1144	// Thus:
1145	//ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
1146	//cout << "PLING" << endl;
1147	/*ring dstRing=rCopy0(srcRing);
1148	for (int ii=0;ii<this->numVars;ii++)
1149	{
1150	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
1151	}*/
1152	ring dstRing=rCopy0(tmpRing);
1153	int length=iv_weight->length();
1154	int A=(int )omAlloc0(length*sizeof(int));
1155	for(int jj=0;jj<length;jj++)
1156	{
1157	A[jj]=(*iv_weight)[jj];
1158	}
1159	dstRing->wvhdl[0]=(int*)A;
1160	rComplete(dstRing);
1161	rChangeCurrRing(dstRing);
1162	//rDelete(tmpRing);
1163	delete iv_weight;
1164	//#ifdef gfan_DEBUG
1165	rWrite(dstRing); cout << endl;
1166	//#endif
1167	ideal dstRing_I;
1168	dstRing_I=idrCopyR(srcRing_HH,srcRing);
1169	//dstRing_I=idrCopyR(inputIdeal,srcRing);
1170	//validOpts<1>=TRUE;
1171	#ifdef gfan_DEBUG
1172	//idShow(dstRing_I);
1173	#endif
1174	BITSET save=test;
1175	test\|=Sy_bit(OPT_REDSB);
1176	test\|=Sy_bit(OPT_REDTAIL);
1177	#ifdef gfan_DEBUG
1178	test\|=Sy_bit(6); //OPT_DEBUG
1179	#endif
1180	ideal tmpI;
1181	//NOTE Any of the two variants of tmpI={idrCopy(),dstRing_I} does the trick
1182	//tmpI = idrCopyR(this->inputIdeal,this->baseRing);
1183	tmpI = dstRing_I;
1184	dstRing_I=kStd(tmpI,NULL,testHomog,NULL);
1185	idNorm(dstRing_I);
1186	//kInterRed(dstRing_I);
1187	idSkipZeroes(dstRing_I);
1188	test=save;
1189	/End of step 3 - reduction/
1190
1191	f->setFlipGB(dstRing_I);//store the flipped GB
1192	f->flipRing=rCopy(dstRing); //store the ring on the other side
1193	//#ifdef gfan_DEBUG
1194	cout << "Flipped GB is UCN " << counter+1 << ":" << endl;
1195	// f->printFlipGB();
1196	this->idDebugPrint(dstRing_I);
1197	cout << endl;
1198	//#endif
1199	rChangeCurrRing(srcRing); //return to the ring we started the computation of flipGB in
1200	// rDelete(dstRing);
1201	}//void flip(ideal gb, facet *f)
1202
1203	/** \brief Compute the remainder of a polynomial by a given ideal
1204	*
1205	* Compute \f$ f^{\mathcal{G}} \f$
1206	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
1207	* However, since we are only interested in the remainder, there is no need to
1208	* compute the factors \f$ a_i \f$
1209	*/
1210	//NOTE: Should be replaced by kNF or kNF2
1211	poly gcone::restOfDiv(poly const &f, ideal const &I)
1212	{
1213	// cout << "Entering restOfDiv" << endl;
1214	poly p=f;
1215	//pWrite(p);
1216	poly r=NULL; //The zero polynomial
1217	int ii;
1218	bool divOccured;
1219
1220	while (p!=NULL)
1221	{
1222	ii=1;
1223	divOccured=FALSE;
1224
1225	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
1226	{
1227	if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
1228	{
1229	poly step1,step2,step3;
1230	//cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
1231	step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
1232	//cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
1233	step2 = ppMult_qq(step1, I->m[ii-1]);
1234	step3 = pSub(pCopy(p), step2);
1235	//p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
1236	//pSetm(p);
1237	pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
1238	p=step3;
1239	//pWrite(p);
1240	divOccured=TRUE;
1241	}
1242	else
1243	{
1244	ii += 1;
1245	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
1246	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
1247	if (divOccured==FALSE)
1248	{
1249	//cout << "TICK 5" << endl;
1250	r=pAdd(pCopy(r),pHead(p));
1251	pSetm(r);
1252	pSort(r);
1253	//cout << "r="; pWrite(r); cout << endl;
1254
1255	if (pLength(p)!=1)
1256	{
1257	p=pSub(pCopy(p),pHead(p)); //Here it may occur that p=0 instead of p=NULL
1258	}
1259	else
1260	{
1261	p=NULL; //Hack to correct this situation
1262	}
1263	//cout << "p="; pWrite(p);
1264	}//if (divOccured==FALSE)
1265	}//while (p!=0)
1266	return r;
1267	}//poly restOfDiv(poly const &f, ideal const &I)
1268
1269	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
1270	*/
1271	//NOTE: use kNF or kNF2 instead of restOfDivision
1272	ideal gcone::ffG(ideal const &H, ideal const &G)
1273	{
1274	// cout << "Entering ffG" << endl;
1275	int size=IDELEMS(H);
1276	ideal res=idInit(size,1);
1277	poly temp1, temp2, temp3; //polys to temporarily store values for pSub
1278	for (int ii=0;ii<size;ii++)
1279	{
1280	// res->m[ii]=restOfDiv(H->m[ii],G);
1281	res->m[ii]=kNF(G, NULL,H->m[ii],0,0);
1282	temp1=H->m[ii];
1283	temp2=res->m[ii];
1284	temp3=pSub(temp1, temp2);
1285	res->m[ii]=temp3;
1286	//res->m[ii]=pSub(temp1,temp2); //buggy
1287	//pSort(res->m[ii]);
1288	//pSetm(res->m[ii]);
1289	//cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
1290	}
1291	return res;
1292	}
1293
1294	/** \brief Compute a Groebner Basis
1295	*
1296	* Computes the Groebner basis and stores the result in gcone::gcBasis
1297	*\param ideal
1298	*\return void
1299	*/
1300	void gcone::getGB(ideal const &inputIdeal)
1301	{
1302	BITSET save=test;
1303	test\|=Sy_bit(OPT_REDSB);
1304	test\|=Sy_bit(OPT_REDTAIL);
1305	ideal gb;
1306	gb=kStd(inputIdeal,NULL,testHomog,NULL);
1307	idNorm(gb);
1308	idSkipZeroes(gb);
1309	this->gcBasis=gb; //write the GB into gcBasis
1310	test=save;
1311	}//void getGB
1312
1313	/** \brief Compute the negative of a given intvec
1314	*/
1315	intvec gcone::ivNeg(const intvec iv)
1316	{
1317	intvec *res = new intvec(iv->length());
1318	res=ivCopy(iv);
1319	res = (int)-1;
1320	return res;
1321	}
1322
1323
1324	/** \brief Compute the dot product of two intvecs
1325	*
1326	*/
1327	int gcone::dotProduct(intvec const &iva, intvec const &ivb)
1328	{
1329	int res=0;
1330	for (int i=0;i<this->numVars;i++)
1331	{
1332	res = res+(iva[i]*ivb[i]);
1333	}
1334	return res;
1335	}//int dotProduct
1336
1337	/** \brief Check whether two intvecs are parallel
1338	*
1339	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1340	*/
1341	bool gcone::isParallel(intvec const &a, intvec const &b)
1342	{
1343	int lhs,rhs;
1344	bool res;
1345	lhs=dotProduct(a,b)*dotProduct(a,b);
1346	rhs=dotProduct(a,a)*dotProduct(b,b);
1347	//cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1348	if (lhs==rhs)
1349	{
1350	res = TRUE;
1351	}
1352	else
1353	{
1354	res = FALSE;
1355	}
1356	return res;
1357	}//bool isParallel
1358
1359	/** \brief Compute an interior point of a given cone
1360	* Result will be written into intvec iv.
1361	* Any rational point is automatically converted into an integer.
1362	*/
1363	void gcone::interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
1364	{
1365	dd_LPPtr lp,lpInt;
1366	dd_ErrorType err=dd_NoError;
1367	dd_LPSolverType solver=dd_DualSimplex;
1368	dd_LPSolutionPtr lpSol=NULL;
1369	dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior
1370	dd_rowindex ddnewpos;
1371	dd_NumberType numb;
1372	//M->representation=dd_Inequality;
1373	//M->objective-dd_LPMin; //Not sure whether this is needed
1374
1375	//NOTE: Make this n-dimensional!
1376	//dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1377
1378	/*NOTE: Leave the following line commented out!
1379	* Otherwise it will cause interior points that are not strictly positive on some examples
1380	*
1381	*/
1382	//dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1383	//if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1384	//cout << "Tick 2" << endl;
1385	//dd_WriteMatrix(stdout,M);
1386
1387	lp=dd_Matrix2LP(M, &err);
1388	if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}
1389	if (lp==NULL){cout << "LP is NULL" << endl;}
1390	#ifdef gfan_DEBUG
1391	// dd_WriteLP(stdout,lp);
1392	#endif
1393
1394	lpInt=dd_MakeLPforInteriorFinding(lp);
1395	if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1396	#ifdef gfan_DEBUG
1397	// dd_WriteLP(stdout,lpInt);
1398	#endif
1399
1400	dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1401	if (err!=dd_NoError)
1402	{
1403	cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1404	dd_WriteErrorMessages(stdout, err);
1405	}
1406
1407	//dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior
1408	if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1409
1410	//lpSol=dd_CopyLPSolution(lpInt);
1411	if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}
1412	#ifdef gfan_DEBUG
1413	cout << "Interior point: ";
1414	for (int ii=1; ii<(lpSol->d)-1;ii++)
1415	{
1416	dd_WriteNumber(stdout,lpSol->sol[ii]);
1417	}
1418	cout << endl;
1419	#endif
1420
1421	//NOTE The following strongly resembles parts of makeInt.
1422	//Maybe merge sometimes
1423	mpz_t kgV; mpz_init(kgV);
1424	mpz_set_str(kgV,"1",10);
1425	mpz_t den; mpz_init(den);
1426	mpz_t tmp; mpz_init(tmp);
1427	mpq_get_den(tmp,lpSol->sol[1]);
1428	for (int ii=1;ii<(lpSol->d)-1;ii++)
1429	{
1430	mpq_get_den(den,lpSol->sol[ii+1]);
1431	mpz_lcm(kgV,tmp,den);
1432	mpz_set(tmp, kgV);
1433	}
1434	mpq_t qkgV;
1435	mpq_init(qkgV);
1436	mpq_set_z(qkgV,kgV);
1437	for (int ii=1;ii<(lpSol->d)-1;ii++)
1438	{
1439	mpq_t product;
1440	mpq_init(product);
1441	mpq_mul(product,qkgV,lpSol->sol[ii]);
1442	iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
1443	mpq_clear(product);
1444	}
1445	#ifdef gfan_DEBUG
1446	// iv.show();
1447	// cout << endl;
1448	#endif
1449	mpq_clear(qkgV);
1450	mpz_clear(tmp);
1451	mpz_clear(den);
1452	mpz_clear(kgV);
1453
1454	dd_FreeLPSolution(lpSol);
1455	dd_FreeLPData(lpInt);
1456	dd_FreeLPData(lp);
1457	set_free(ddlinset);
1458	set_free(ddredrows);
1459
1460	}//void interiorPoint(dd_MatrixPtr const &M)
1461
1462	/** \brief Copy a ring and add a weighted ordering in first place
1463	*
1464	*/
1465	ring gcone::rCopyAndAddWeight(ring const &r, intvec const *ivw)
1466	{
1467	ring res=rCopy0(r);
1468	int jj;
1469
1470	omFree(res->order);
1471	res->order =(int )omAlloc0(4sizeof(int));
1472	omFree(res->block0);
1473	res->block0=(int )omAlloc0(4sizeof(int));
1474	omFree(res->block1);
1475	res->block1=(int )omAlloc0(4sizeof(int));
1476	omfree(res->wvhdl);
1477	res->wvhdl =(int *)omAlloc0(4sizeof(int**));
1478
1479	res->order[0]=ringorder_a;
1480	res->block0[0]=1;
1481	res->block1[0]=res->N;;
1482	res->order[1]=ringorder_dp; //basically useless, since that should never be used
1483	res->block0[1]=1;
1484	res->block1[1]=res->N;;
1485	res->order[2]=ringorder_C;
1486
1487	int length=ivw->length();
1488	int A=(int )omAlloc0(length*sizeof(int));
1489	for (jj=0;jj<length;jj++)
1490	{
1491	A[jj]=(*ivw)[jj];
1492	}
1493	res->wvhdl[0]=(int *)A;
1494	res->block1[0]=length;
1495
1496	rComplete(res);
1497	return res;
1498	}//rCopyAndAdd
1499
1500	ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1501	{
1502	ring res=rCopy0(currRing);
1503	rComplete(res);
1504	rSetWeightVec(res,(int64*)ivw);
1505	//rChangeCurrRing(rnew);
1506	return res;
1507	}
1508
1509	/** \brief Checks whether a given facet is a search facet
1510	* Determines whether a given facet of a cone is the search facet of a neighbouring cone
1511	* This is done in the following way:
1512	* We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1513	* that is first crossed during the generic walk.
1514	* We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1515	* If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1516	*/
1517	bool gcone::isSearchFacet(gcone &gcTmp, facet *testfacet)
1518	{
1519	ring actRing=currRing;
1520	facet facetPtr=(facet)gcTmp.facetPtr;
1521	facet fMin=new facet(facetPtr); //Pointer to the "minimal" facet
1522	//facet *fMin = new facet(tmpcone.facetPtr);
1523	//fMin=tmpcone.facetPtr; //Initialise to first facet of tmpcone
1524	facet *fAct; //Ptr to alpha_i
1525	facet *fCmp; //Ptr to alpha_j
1526	fAct = fMin;
1527	fCmp = fMin->next;
1528
1529	rChangeCurrRing(this->rootRing); //because we compare the monomials in the rootring
1530	poly p=pInit();
1531	poly q=pInit();
1532	intvec *alpha_i = new intvec(this->numVars);
1533	intvec *alpha_j = new intvec(this->numVars);
1534	intvec *sigma = new intvec(this->numVars);
1535	sigma=gcTmp.getIntPoint();
1536
1537	int u=(int )omAlloc((this->numVars+1)*sizeof(int));
1538	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1539	u[0]=0; v[0]=0;
1540	int weight1,weight2;
1541	while(fAct->next->next!=NULL) //NOTE this is ugly. Can it be done without fCmp?
1542	{
1543	/* Get alpha_i and alpha_{i+1} */
1544	alpha_i=fAct->getFacetNormal();
1545	alpha_j=fCmp->getFacetNormal();
1546	#ifdef gfan_DEBUG
1547	alpha_i->show();
1548	alpha_j->show();
1549	#endif
1550	/Compute the dot product of sigma and alpha_{i,j}/
1551	weight1=dotProduct(sigma,alpha_i);
1552	weight2=dotProduct(sigma,alpha_j);
1553	#ifdef gfan_DEBUG
1554	cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1555	#endif
1556	/Adjust alpha_i and alpha_i+1 accordingly/
1557	for(int ii=1;ii<=this->numVars;ii++)
1558	{
1559	u[ii]=weight1(alpha_i)[ii-1];
1560	v[ii]=weight2(alpha_j)[ii-1];
1561	}
1562
1563	/Now p_weight and q_weight need to be compared as exponent vectors/
1564	pSetCoeff0(p,nInit(1));
1565	pSetCoeff0(q,nInit(1));
1566	pSetExpV(p,u);
1567	pSetm(p);
1568	pSetExpV(q,v);
1569	pSetm(q);
1570	#ifdef gfan_DEBUG
1571	pWrite(p);pWrite(q);
1572	#endif
1573	/*We want to check whether x^p < x^q
1574	=> want to check for return value 1 */
1575	if (pLmCmp(p,q)==1) //i.e. x^q is smaller
1576	{
1577	fMin=fCmp;
1578	fAct=fMin;
1579	fCmp=fCmp->next;
1580	}
1581	else
1582	{
1583	//fAct=fAct->next;
1584	if(fCmp->next!=NULL)
1585	{
1586	fCmp=fCmp->next;
1587	}
1588	else
1589	{
1590	fAct=fAct->next;
1591	}
1592	}
1593	//fAct=fAct->next;
1594	}//while(fAct.facetPtr->next!=NULL)
1595	delete alpha_i,alpha_j,sigma;
1596
1597	/If testfacet was minimal then fMin should still point there /
1598
1599	//if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))
1600	#ifdef gfan_DEBUG
1601	cout << "Checking for parallelity" << endl <<" fMin is";
1602	fMin->printNormal();
1603	cout << "testfacet is ";
1604	testfacet->printNormal();
1605	cout << endl;
1606	#endif
1607	if (fMin==gcTmp.facetPtr)
1608	//if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
1609	//if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
1610	{
1611	cout << "Parallel" << endl;
1612	rChangeCurrRing(actRing);
1613	//delete alpha_min, test;
1614	return TRUE;
1615	}
1616	else
1617	{
1618	cout << "Not parallel" << endl;
1619	rChangeCurrRing(actRing);
1620	//delete alpha_min, test;
1621	return FALSE;
1622	}
1623	}//bool isSearchFacet
1624
1625	/** \brief Check for equality of two intvecs
1626	*/
1627	bool gcone::areEqual(intvec const &a, intvec const &b)
1628	{
1629	bool res=TRUE;
1630	for(int ii=0;ii<this->numVars;ii++)
1631	{
1632	if(a[ii]!=b[ii])
1633	{
1634	res=FALSE;
1635	break;
1636	}
1637	}
1638	return res;
1639	}
1640
1641	/** \brief The reverse search algorithm
1642	*/
1643	void gcone::reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1644	{
1645	facet *fAct=new facet();
1646	fAct = gcAct->facetPtr;
1647
1648	while(fAct!=NULL)
1649	{
1650	cout << "======================"<< endl;
1651	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1652	gcone gcTmp = new gcone(gcAct);
1653	//idShow(gcTmp->gcBasis);
1654	gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1655	#ifdef gfan_DEBUG
1656	facet *f = new facet();
1657	f=gcTmp->facetPtr;
1658	while(f!=NULL)
1659	{
1660	f->printNormal();
1661	f=f->next;
1662	}
1663	#endif
1664	gcTmp->showIntPoint();
1665	/recursive part goes gere/
1666	if (isSearchFacet(gcTmp,(facet)gcAct->facetPtr))
1667	{
1668	gcAct->next=gcTmp;
1669	cout << "PING"<< endl;
1670	reverseSearch(gcTmp);
1671	}
1672	else
1673	{
1674	delete gcTmp;
1675	/NOTE remove fAct from linked list. It's no longer needed/
1676	}
1677	/recursion ends/
1678	fAct = fAct->next;
1679	}//while(fAct->next!=NULL)
1680	}//reverseSearch
1681
1682	/** \brief The new method of Markwig and Jensen
1683	* Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1684	* In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1685	* Keep a list of facets as a linked list containing an intvec and an integer matrix.
1686	* Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1687	* soon as we compute this facet again. Comparison of facets is done by...
1688	*/
1689	void gcone::noRevS(gcone &gcRoot, bool usingIntPoint)
1690	{
1691	facet *SearchListRoot = new facet(); //The list containing ALL facets we come across
1692	facet *SearchListAct;
1693	SearchListAct = NULL;
1694	//SearchListAct = SearchListRoot;
1695
1696	gcone *gcAct;
1697	gcAct = &gcRoot;
1698	gcone *gcPtr; //Pointer to end of linked list of cones
1699	gcPtr = &gcRoot;
1700	gcone *gcNext; //Pointer to next cone to be visited
1701	gcNext = &gcRoot;
1702	gcone *gcHead;
1703	gcHead = &gcRoot;
1704
1705	facet *fAct;
1706	fAct = gcAct->facetPtr;
1707
1708	ring rAct;
1709	rAct = currRing;
1710
1711	int UCNcounter=gcAct->getUCN();
1712
1713	/* Use pure SLA or writeCone2File */
1714	// enum heuristic {
1715	// ram,
1716	// hdd
1717	// };
1718	// heuristic h;
1719	// h=hdd;
1720
1721	#ifdef gfan_DEBUG
1722	cout << "NoRevs" << endl;
1723	cout << "Facets are:" << endl;
1724	gcAct->showFacets();
1725	#endif
1726
1727	gcAct->getCodim2Normals(*gcAct);
1728
1729	//Compute unique representation of codim-2-facets
1730	gcAct->normalize();
1731
1732	/* Make a copy of the facet list of first cone
1733	Since the operations getCodim2Normals and normalize affect the facets
1734	we must not memcpy them before these ops!
1735	*/
1736	intvec *fNormal;// = new intvec(this->numVars);
1737	intvec *f2Normal;// = new intvec(this->numVars);
1738	facet *codim2Act; codim2Act = NULL;
1739	facet *sl2Root; //sl2Root = new facet(2);
1740	facet *sl2Act; //sl2Act = sl2Root;
1741
1742	for(int ii=0;ii<this->numFacets;ii++)
1743	{
1744	//only copy flippable facets! NOTE: We do not need flipRing or any such information.
1745	if(fAct->isFlippable==TRUE)
1746	{
1747	fNormal = fAct->getFacetNormal();
1748	if( ii==0 \|\| (ii>0 && SearchListAct==NULL) ) //1st facet may be non-flippable
1749	{
1750	SearchListAct = SearchListRoot;
1751	//memcpy(SearchListAct,fAct,sizeof(facet));
1752	}
1753	else
1754	{
1755	SearchListAct->next = new facet();
1756	SearchListAct = SearchListAct->next;
1757	//memcpy(SearchListAct,fAct,sizeof(facet));
1758	}
1759	SearchListAct->setFacetNormal(fNormal);
1760	SearchListAct->setUCN(this->getUCN());
1761	SearchListAct->numCodim2Facets=fAct->numCodim2Facets;
1762	SearchListAct->isFlippable=TRUE;
1763	//Copy codim2-facets
1764	codim2Act=fAct->codim2Ptr;
1765	SearchListAct->codim2Ptr = new facet(2);
1766	sl2Root = SearchListAct->codim2Ptr;
1767	sl2Act = sl2Root;
1768	//while(codim2Act!=NULL)
1769	for(int jj=0;jj<fAct->numCodim2Facets;jj++)
1770	{
1771	f2Normal = codim2Act->getFacetNormal();
1772	if(jj==0)
1773	{
1774	sl2Act = sl2Root;
1775	sl2Act->setFacetNormal(f2Normal);
1776	}
1777	else
1778	{
1779	sl2Act->next = new facet(2);
1780	sl2Act = sl2Act->next;
1781	sl2Act->setFacetNormal(f2Normal);
1782	}
1783	codim2Act = codim2Act->next;
1784	}
1785	fAct = fAct->next;
1786	}//if(fAct->isFlippable==TRUE)
1787	else {fAct = fAct->next;}
1788	}//End of copying facets into SLA
1789
1790	SearchListAct = SearchListRoot; //Set to beginning of list
1791	/Make SearchList doubly linked/
1792	while(SearchListAct!=NULL)
1793	{
1794	if(SearchListAct->next!=NULL)
1795	{
1796	SearchListAct->next->prev = SearchListAct;
1797	}
1798	SearchListAct = SearchListAct->next;
1799	}
1800	SearchListAct = SearchListRoot; //Set to beginning of List
1801
1802	fAct = gcAct->facetPtr;
1803	//NOTE Disabled until code works as expected
1804	//gcAct->writeConeToFile(*gcAct);
1805
1806	/End of initialisation/
1807	fAct = SearchListAct;
1808	/*2nd step
1809	Choose a facet from fListPtr, flip it and forget the previous cone
1810	We always choose the first facet from fListPtr as facet to be flipped
1811	*/
1812	while((SearchListAct!=NULL))// && counter<10)
1813	{//NOTE See to it that the cone is only changed after ALL facets have been flipped!
1814	fAct = SearchListAct;
1815
1816	while(fAct!=NULL)
1817	{ //Since SLA should only contain flippables there should be no need to check for that
1818	gcAct->flip(gcAct->gcBasis,fAct);
1819	ring rTmp=rCopy(fAct->flipRing);
1820	rComplete(rTmp);
1821	rChangeCurrRing(rTmp);
1822	gcone gcTmp = new gcone::gcone(gcAct,*fAct);
1823	/* Now gcTmp->gcBasis and gcTmp->baseRing are set from fAct->flipGB and fAct->flipRing.
1824	* Since we come at most once across a given facet from gcAct->facetPtr we can delete them.
1825	* NOTE: Can this cause trouble with the destructor?
1826	* Answer: Yes it bloody well does!
1827	* However I'd like to delete this data here, since if we have a cone with many many facets it
1828	* should pay to get rid of the flipGB as soon as possible.
1829	* Destructor must be equipped with necessary checks.
1830	*/
1831	// idDelete((ideal *)&fAct->flipGB);
1832	// rDelete(fAct->flipRing);
1833	gcTmp->getConeNormals(gcTmp->gcBasis, FALSE);
1834	gcTmp->getCodim2Normals(*gcTmp);
1835	gcTmp->normalize();
1836	#ifdef gfan_DEBUG
1837	gcTmp->showFacets(1);
1838	#endif
1839	if(gfanHeuristic==1)
1840	{
1841	gcTmp->writeConeToFile(*gcTmp);
1842	}
1843	/add facets to SLA here/
1844	SearchListRoot=gcTmp->enqueueNewFacets(*SearchListRoot);
1845	#ifdef gfan_DEBUG
1846	if(SearchListRoot!=NULL)
1847	gcTmp->showSLA(*SearchListRoot);
1848	#endif
1849	rChangeCurrRing(gcAct->baseRing);
1850	//rDelete(rTmp);
1851	//doubly linked for easier removal
1852	gcTmp->prev = gcPtr;
1853	gcPtr->next=gcTmp;
1854	gcPtr=gcPtr->next;
1855	if(fAct->getUCN() == fAct->next->getUCN())
1856	{
1857	fAct=fAct->next;
1858	}
1859	else
1860	break;
1861	}//while( ( (fAct->next!=NULL) && (fAct->getUCN()==fAct->next->getUCN() ) ) );
1862	//NOTE Now all neighbouring cones of gcAct have been computed, so we may delete gcAct
1863	gcone *deleteMarker;
1864	deleteMarker=gcAct;
1865	// delete deleteMarker;
1866	// deleteMarker=NULL;
1867	//Search for cone with smallest UCN
1868	gcNext = gcHead;
1869	while(gcNext!=NULL && SearchListRoot!=NULL && gcNext!=(gcone * const)0xfbfbfbfbfbfbfbfb && gcNext!=(gcone * const)0xfbfbfbfb)
1870	{
1871	/*NOTE if gcNext->getUCN is smaller than SearchListRoot->getUCN it follows, that the cone
1872	gcNext will not be used in any more computation. So -> delete
1873	*/
1874	if (gcNext->getUCN() < SearchListRoot->getUCN() )
1875	{
1876	//idDelete((ideal *)&gcNext->gcBasis); //at least drop huge grÃ¶bner bases
1877	if(gcNext == gcHead)
1878	{
1879	deleteMarker = gcHead;
1880	gcHead = gcNext->next;
1881	//gcNext->next->prev = NULL;
1882	}
1883	else
1884	{
1885	deleteMarker = gcNext;
1886	gcNext->prev->next = gcNext->next;
1887	gcNext->next->prev = gcNext->prev;
1888	}
1889	// gcNext = gcNext->next;
1890	// delete deleteMarker;
1891	// deleteMarker = NULL;
1892	}
1893	else if( gcNext->getUCN() == SearchListRoot->getUCN() )
1894	{//NOTE: Implement heuristic to be used!
1895	if(gfanHeuristic==0)
1896	{
1897	gcAct = gcNext;
1898	rAct=rCopy(gcAct->baseRing);
1899	rComplete(rAct);
1900	rChangeCurrRing(rAct);
1901	break;
1902	}
1903	else if(gfanHeuristic==1)
1904	{
1905	//Read st00f from file
1906	gcAct = gcNext;
1907	//implant the GB into gcAct
1908	readConeFromFile(gcNext->getUCN(), gcAct);
1909	rAct=rCopy(gcAct->baseRing);
1910	rComplete(rAct);
1911	rChangeCurrRing(rAct);
1912	break;
1913	}
1914	}
1915	gcNext = gcNext->next;
1916	// if(deleteMarker!=NULL && deleteMarker!=(gcone *)0xfbfbfbfbfbfbfbfb)
1917	if(this->getUCN()!=1) //workaround to avoid unpredictable behaviour towards end of program
1918	delete deleteMarker;
1919	// deleteMarker = NULL;
1920	// gcNext = gcNext->next;
1921	}
1922	UCNcounter++;
1923	//SearchListAct = SearchListAct->next;
1924	//SearchListAct = fAct->next;
1925	SearchListAct = SearchListRoot;
1926	}
1927	cout << endl << "Found " << counter << " cones - terminating" << endl;
1928
1929	//NOTE Hm, comment in and get a crash for free...
1930	//dd_free_global_constants();
1931	//gc.writeConeToFile(gc);
1932
1933
1934	//fAct = fListPtr;
1935	//gcone *gcTmp = new gcone(gc); //copy
1936	//gcTmp->flip(gcTmp->gcBasis,fAct);
1937	//NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
1938
1939	}//void noRevS(gcone &gc)
1940
1941
1942	/** \brief Make a set of rational vectors into integers
1943	*
1944	* We compute the lcm of the denominators and multiply with this to get integer values.
1945	* \param dd_MatrixPtr,intvec
1946	*/
1947	void gcone::makeInt(dd_MatrixPtr const &M, int const line, intvec &n)
1948	{
1949	mpz_t denom[this->numVars];
1950	for(int ii=0;ii<this->numVars;ii++)
1951	{
1952	mpz_init_set_str(denom[ii],"0",10);
1953	}
1954
1955	mpz_t kgV,tmp;
1956	mpz_init(kgV);
1957	mpz_init(tmp);
1958
1959	for (int ii=0;ii<(M->colsize)-1;ii++)
1960	{
1961	mpz_t z;
1962	mpz_init(z);
1963	mpq_get_den(z,M->matrix[line-1][ii+1]);
1964	mpz_set( denom[ii], z);
1965	mpz_clear(z);
1966	}
1967
1968	/Compute lcm of the denominators/
1969	mpz_set(tmp,denom[0]);
1970	for (int ii=0;ii<(M->colsize)-1;ii++)
1971	{
1972	mpz_lcm(kgV,tmp,denom[ii]);
1973	mpz_set(tmp,kgV);
1974	}
1975
1976	/Multiply the nominators by kgV/
1977	mpq_t qkgV,res;
1978	mpq_init(qkgV);
1979	mpq_set_str(qkgV,"1",10);
1980	mpq_canonicalize(qkgV);
1981
1982	mpq_init(res);
1983	mpq_set_str(res,"1",10);
1984	mpq_canonicalize(res);
1985
1986	mpq_set_num(qkgV,kgV);
1987
1988	// mpq_canonicalize(qkgV);
1989	for (int ii=0;ii<(M->colsize)-1;ii++)
1990	{
1991	mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
1992	n[ii]=(int)mpz_get_d(mpq_numref(res));
1993	}
1994	//mpz_clear(denom[this->numVars]);
1995	mpz_clear(kgV);
1996	mpq_clear(qkgV); mpq_clear(res);
1997
1998	}
1999	/**
2000	* For all codim-2-facets we compute the gcd of the components of the facet normal and
2001	* divide it out. Thus we get a normalized representation of each
2002	* (codim-2)-facet normal, i.e. a primitive vector
2003	*/
2004	void gcone::normalize()
2005	{
2006	int ggT=1;
2007	facet *fAct;
2008	facet *codim2Act;
2009	fAct = this->facetPtr;
2010	codim2Act = fAct->codim2Ptr;
2011	intvec *n = new intvec(this->numVars);
2012
2013	//while(codim2Act->next!=NULL)
2014	while(fAct!=NULL)
2015	{
2016	while(codim2Act!=NULL)
2017	{
2018	n=codim2Act->getFacetNormal();
2019	for(int ii=0;ii<this->numVars;ii++)
2020	{
2021	ggT = intgcd(ggT,(*n)[ii]);
2022	}
2023	for(int ii=0;ii<this->numVars;ii++)
2024	{
2025	(n)[ii] = ((n)[ii])/ggT;
2026	}
2027	codim2Act->setFacetNormal(n);
2028	codim2Act = codim2Act->next;
2029	}
2030	fAct = fAct->next;
2031	}
2032	delete n;
2033
2034	}
2035	/** \brief Enqueue new facets into the searchlist
2036	* The searchlist (SLA for short) is implemented as a FIFO
2037	* Checks are done as follows:
2038	* 1) Check facet fAct against all facets in SLA for parallelity. If it is not parallel to any one add it.
2039	* 2) If it is parallel compare the codim2 facets. If they coincide the facets are equal and we delete the
2040	* facet from SLA and do not add fAct.
2041	* It may very well happen, that SLA=\f$ \emptyset \f$ but we do not have checked all fActs yet. In this case we
2042	* can be sure, that none of the facets that are still there already were in SLA once, so we just dump them into SLA.
2043	* Takes ptr to search list root
2044	* Returns a pointer to new first element of Searchlist
2045	*/
2046	//void enqueueNewFacets(facet &f)
2047	facet * gcone::enqueueNewFacets(facet &f)
2048	{
2049	facet *slHead;
2050	slHead = &f;
2051	facet *slAct; //called with f=SearchListRoot
2052	slAct = &f;
2053	facet *slEnd; //Pointer to end of SLA
2054	slEnd = &f;
2055	facet *slEndStatic; //marks the end before a new facet is added
2056	facet *fAct;
2057	fAct = this->facetPtr;
2058	facet *codim2Act;
2059	codim2Act = this->facetPtr->codim2Ptr;
2060	facet *sl2Act;
2061	sl2Act = f.codim2Ptr;
2062	/** Pointer to a facet that will be deleted */
2063	facet *deleteMarker;
2064	deleteMarker = NULL;
2065
2066	/** Flag to indicate a facet that should be added to SLA*/
2067	bool doNotAdd=FALSE;
2068	/** \brief Flag to mark a facet that might be added
2069	* The following case may occur:
2070	* A facet fAct is found to be parallel but not equal to the current facet slAct from SLA.
2071	* This does however not mean that there does not exist a facet behind the current slAct that is actually equal
2072	* to fAct. In this case we set the boolean flag maybe to TRUE. If we encounter an equality lateron, it is reset to
2073	* FALSE and the according slAct is deleted.
2074	* If slAct->next==NULL AND maybe==TRUE we know, that fAct must be added
2075	*/
2076	volatile bool maybe=FALSE;
2077	/**A facet was removed, lengthOfSearchlist-- thus we must not rely on
2078	* if(notParallelCtr==lengthOfSearchList) but rather
2079	* if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE)
2080	*/
2081	volatile bool removalOccured=FALSE;
2082	int ctr=0; //encountered equalities in SLA
2083	int notParallelCtr=0;
2084	int lengthOfSearchList=1;
2085	while(slEnd->next!=NULL)
2086	{
2087	slEnd=slEnd->next;
2088	lengthOfSearchList++;
2089	}
2090	slEndStatic = slEnd;
2091	/1st step: compare facetNormals/
2092	intvec *fNormal=NULL; //= new intvec(this->numVars);
2093	intvec *f2Normal=NULL; //= new intvec(this->numVars);
2094	intvec *slNormal=NULL; //= new intvec(this->numVars);
2095	intvec *sl2Normal=NULL; //= new intvec(this->numVars);
2096
2097	while(fAct!=NULL)
2098	{
2099	if(fAct->isFlippable==TRUE)
2100	{
2101	maybe=FALSE;
2102	doNotAdd=TRUE;
2103	fNormal=fAct->getFacetNormal();
2104	slAct = slHead;
2105	notParallelCtr=0;
2106	// delete deleteMarker;
2107	// deleteMarker=NULL;
2108	/*If slAct==NULL and fAct!=NULL
2109	we just copy all remaining facets into SLA*/
2110	if(slAct==NULL)
2111	{
2112	facet *fCopy;
2113	fCopy = fAct;
2114	while(fCopy!=NULL)
2115	{
2116	if(slAct==NULL)
2117	{
2118	slAct = new facet();
2119	slHead = slAct;
2120	}
2121	else
2122	{
2123	slAct->next = new facet();
2124	slAct = slAct->next;
2125	}
2126	slAct->setFacetNormal(fAct->getFacetNormal());
2127	slAct->setUCN(fAct->getUCN());
2128	slAct->isFlippable=fAct->isFlippable;
2129	facet *f2Copy;
2130	f2Copy = fCopy->codim2Ptr;
2131	sl2Act = slAct->codim2Ptr;
2132	while(f2Copy!=NULL)
2133	{
2134	if(sl2Act==NULL)
2135	{
2136	sl2Act = new facet(2);
2137	slAct->codim2Ptr = sl2Act;
2138	}
2139	else
2140	{
2141	facet *marker;
2142	marker = sl2Act;
2143	sl2Act->next = new facet(2);
2144	sl2Act = sl2Act->next;
2145	sl2Act->prev = marker;
2146	}
2147	sl2Act->setFacetNormal(f2Copy->getFacetNormal());
2148	sl2Act->setUCN(f2Copy->getUCN());
2149	f2Copy = f2Copy->next;
2150	}
2151	fCopy = fCopy->next;
2152	}
2153	break;
2154	}
2155	/End of dumping into SLA/
2156	while(slAct!=NULL)
2157	//while(slAct!=slEndStatic->next)
2158	{
2159	// if(deleteMarker!=NULL)
2160	// {
2161	// delete deleteMarker;
2162	// deleteMarker=NULL;
2163	// }
2164	removalOccured=FALSE;
2165	slNormal = slAct->getFacetNormal();
2166	#ifdef gfan_DEBUG
2167	cout << "Checking facet (";
2168	fNormal->show(1,1);
2169	cout << ") against (";
2170	slNormal->show(1,1);
2171	cout << ")" << endl;
2172	#endif
2173	if(!isParallel(fNormal, slNormal))
2174	{
2175	notParallelCtr++;
2176	// slAct = slAct->next;
2177	}
2178	else //fN and slN are parallel
2179	{
2180	//We check whether the codim-2-facets coincide
2181	codim2Act = fAct->codim2Ptr;
2182	ctr=0;
2183	while(codim2Act!=NULL)
2184	{
2185	f2Normal = codim2Act->getFacetNormal();
2186	//sl2Act = f.codim2Ptr;
2187	sl2Act = slAct->codim2Ptr;
2188	while(sl2Act!=NULL)
2189	{
2190	sl2Normal = sl2Act->getFacetNormal();
2191	if(areEqual(f2Normal, sl2Normal))
2192	ctr++;
2193	sl2Act = sl2Act->next;
2194	}//while(sl2Act!=NULL)
2195	codim2Act = codim2Act->next;
2196	}//while(codim2Act!=NULL)
2197	if(ctr==fAct->numCodim2Facets) //facets ARE equal
2198	{
2199	#ifdef gfan_DEBUG
2200	cout << "Removing facet (";
2201	slNormal->show(1,0);
2202	cout << ") from SLA:" << endl;
2203	fAct->fDebugPrint();
2204	slAct->fDebugPrint();
2205	#endif
2206	removalOccured=TRUE;
2207	slAct->isFlippable=FALSE;
2208	doNotAdd=TRUE;
2209	maybe=FALSE;
2210	if(slAct==slHead) //We want to delete the first element of SearchList
2211	{
2212	deleteMarker = slHead;
2213	slHead = slAct->next;
2214	if(slHead!=NULL)
2215	slHead->prev = NULL;
2216	//set a bool flag to mark slAct as to be deleted
2217	}//NOTE find a way to delete without affecting slAct = slAct->next
2218	else if(slAct==slEndStatic)
2219	{
2220	deleteMarker = slAct;
2221	if(slEndStatic->next==NULL)
2222	{
2223	slEndStatic = slEndStatic->prev;
2224	slEndStatic->next = NULL;
2225	slEnd = slEndStatic;
2226	}
2227	else //we already added a facet after slEndStatic
2228	{
2229	slEndStatic->prev->next = slEndStatic->next;
2230	slEndStatic->next->prev = slEndStatic->prev;
2231	slEndStatic = slEndStatic->prev;
2232	//slEnd = slEndStatic;
2233	}
2234	}
2235	else
2236	{
2237	deleteMarker = slAct;
2238	slAct->prev->next = slAct->next;
2239	slAct->next->prev = slAct->prev;
2240	}
2241
2242	//update lengthOfSearchList
2243	lengthOfSearchList--;
2244	//delete slAct;
2245	//slAct=NULL;
2246	//slAct = slAct->next; //not needed, since facets are equal
2247	//delete deleteMarker;
2248	//deleteMarker=NULL;
2249	//fAct = fAct->next;
2250	break;
2251	}//if(ctr==fAct->numCodim2Facets)
2252	else //facets are parallel but NOT equal. But this does not imply that there
2253	//is no other facet later in SLA that might be equal.
2254	{
2255	maybe=TRUE;
2256	// if(slAct->next==NULL && maybe==TRUE)
2257	// {
2258	// doNotAdd=FALSE;
2259	// slAct = slAct->next;
2260	// break;
2261	// }
2262	// else
2263	// slAct=slAct->next;
2264	}
2265	//slAct = slAct->next;
2266	//delete deleteMarker;
2267	}//if(!isParallel(fNormal, slNormal))
2268	if( (slAct->next==NULL) && (maybe==TRUE) )
2269	{
2270	doNotAdd=FALSE;
2271	break;
2272	}
2273	slAct = slAct->next;
2274	/* NOTE The following lines must not be here but rather called inside the if-block above.
2275	Otherwise results get crappy. Unfortunately there are two slAct=slAct->next calls now,
2276	(not nice!) but since they are in seperate branches of the if-statement there should not
2277	be a way it gets called twice thus ommiting one facet:
2278	slAct = slAct->next;*/
2279	//delete deleteMarker;
2280	//deleteMarker=NULL;
2281	//if slAct was marked as to be deleted, delete it here!
2282	}//while(slAct!=NULL)
2283	if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) \|\| (doNotAdd==FALSE) )
2284	//if( (notParallelCtr==lengthOfSearchList ) \|\| doNotAdd==FALSE )
2285	{
2286	#ifdef gfan_DEBUG
2287	cout << "Adding facet (";
2288	fNormal->show(1,0);
2289	cout << ") to SLA " << endl;
2290	#endif
2291	//Add fAct to SLA
2292	facet *marker;
2293	marker = slEnd;
2294	facet *f2Act;
2295	f2Act = fAct->codim2Ptr;
2296
2297	slEnd->next = new facet();
2298	slEnd = slEnd->next;
2299	facet *slEndCodim2Root;
2300	facet *slEndCodim2Act;
2301	slEndCodim2Root = slEnd->codim2Ptr;
2302	slEndCodim2Act = slEndCodim2Root;
2303
2304	slEnd->setUCN(this->getUCN());
2305	slEnd->isFlippable = TRUE;
2306	slEnd->setFacetNormal(fNormal);
2307	slEnd->prev = marker;
2308	//Copy codim2-facets
2309	intvec *f2Normal;// = new intvec(this->numVars);
2310	while(f2Act!=NULL)
2311	{
2312	f2Normal=f2Act->getFacetNormal();
2313	if(slEndCodim2Root==NULL)
2314	{
2315	slEndCodim2Root = new facet(2);
2316	slEnd->codim2Ptr = slEndCodim2Root;
2317	slEndCodim2Root->setFacetNormal(f2Normal);
2318	slEndCodim2Act = slEndCodim2Root;
2319	}
2320	else
2321	{
2322	slEndCodim2Act->next = new facet(2);
2323	slEndCodim2Act = slEndCodim2Act->next;
2324	slEndCodim2Act->setFacetNormal(f2Normal);
2325	}
2326	f2Act = f2Act->next;
2327	}
2328	lengthOfSearchList++;
2329	}//if( (notParallelCtr==lengthOfSearchList && removalOccured==FALSE) \|\|
2330	fAct = fAct->next;
2331	}//if(fAct->isFlippable==TRUE)
2332	else
2333	{
2334	fAct = fAct->next;
2335	}
2336	}//while(fAct!=NULL)
2337	return slHead;
2338	// delete sl2Normal;
2339	// delete slNormal;
2340	// delete f2Normal;
2341	// delete fNormal;
2342	}//addC2N
2343
2344	/** \brief Compute the gcd of two ints
2345	*/
2346	int gcone::intgcd(int a, int b)
2347	{
2348	int r, p=a, q=b;
2349	if(p < 0)
2350	p = -p;
2351	if(q < 0)
2352	q = -q;
2353	while(q != 0)
2354	{
2355	r = p % q;
2356	p = q;
2357	q = r;
2358	}
2359	return p;
2360	}
2361
2362	/** \brief Construct a dd_MatrixPtr from a cone's list of facets
2363	*
2364	*/
2365	dd_MatrixPtr gcone::facets2Matrix(gcone const &gc)
2366	{
2367	facet *fAct;
2368	fAct = gc.facetPtr;
2369
2370	dd_MatrixPtr M;
2371	dd_rowrange ddrows;
2372	dd_colrange ddcols;
2373	ddcols=(this->numVars)+1;
2374	ddrows=this->numFacets;
2375	dd_NumberType numb = dd_Integer;
2376	M=dd_CreateMatrix(ddrows,ddcols);
2377
2378	int jj=0;
2379
2380	while(fAct!=NULL)
2381	{
2382	intvec *comp;
2383	comp = fAct->getFacetNormal();
2384	for(int ii=0;ii<this->numVars;ii++)
2385	{
2386	dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
2387	}
2388	jj++;
2389	fAct=fAct->next;
2390	}
2391
2392	return M;
2393	}
2394
2395	/** \brief Write information about a cone into a file on disk
2396	*
2397	* This methods writes the information needed for the "second" method into a file.
2398	* The file's is divided in sections containing information on
2399	* 1) the ring
2400	* 2) the cone's Groebner Basis
2401	* 3) the cone's facets
2402	* Each line contains exactly one date
2403	* Each section starts with its name in CAPITALS
2404	*/
2405	void gcone::writeConeToFile(gcone const &gc, bool usingIntPoints)
2406	{
2407	int UCN=gc.UCN;
2408	stringstream ss;
2409	ss << UCN;
2410	string UCNstr = ss.str();
2411
2412	string prefix="/tmp/cone";
2413	string suffix=".sg";
2414	string filename;
2415	filename.append(prefix);
2416	filename.append(UCNstr);
2417	filename.append(suffix);
2418
2419	ofstream gcOutputFile(filename.c_str());
2420	facet *fAct;
2421	fAct = gc.facetPtr;
2422	facet *f2Act;
2423	f2Act=fAct->codim2Ptr;
2424
2425	char *ringString = rString(currRing);
2426
2427	if (!gcOutputFile)
2428	{
2429	cout << "Error opening file for writing in writeConeToFile" << endl;
2430	}
2431	else
2432	{
2433	gcOutputFile << "UCN" << endl;
2434	gcOutputFile << this->UCN << endl;
2435	gcOutputFile << "RING" << endl;
2436	gcOutputFile << ringString << endl;
2437	gcOutputFile << "GCBASISLENGTH" << endl;
2438	gcOutputFile << IDELEMS(gc.gcBasis) << endl;
2439	//Write this->gcBasis into file
2440	gcOutputFile << "GCBASIS" << endl;
2441	for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
2442	{
2443	gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
2444	}
2445
2446	gcOutputFile << "FACETS" << endl;
2447
2448	while(fAct!=NULL)
2449	{
2450	intvec *iv;
2451	intvec *iv2;
2452	iv=fAct->getFacetNormal();
2453	f2Act=fAct->codim2Ptr;
2454	for (int ii=0;ii<iv->length();ii++)
2455	{
2456	if (ii<iv->length()-1)
2457	{
2458	gcOutputFile << (*iv)[ii] << ",";
2459	}
2460	else
2461	{
2462	gcOutputFile << (*iv)[ii] << " ";
2463	}
2464	}
2465	while(f2Act!=NULL)
2466	{
2467	iv2=f2Act->getFacetNormal();
2468	for(int jj=0;jj<iv2->length();jj++)
2469	{
2470	if (jj<iv2->length()-1)
2471	{
2472	gcOutputFile << (*iv2)[jj] << ",";
2473	}
2474	else
2475	{
2476	gcOutputFile << (*iv2)[jj] << " ";
2477	}
2478	}
2479	f2Act = f2Act->next;
2480	}
2481	gcOutputFile << endl;
2482	fAct=fAct->next;
2483	}
2484	gcOutputFile.close();
2485	}
2486
2487	}//writeConeToFile(gcone const &gc)
2488
2489	/** \brief Reads a cone from a file identified by its number
2490	*/
2491	void gcone::readConeFromFile(int UCN, gcone *gc)
2492	{
2493	//int UCN=gc.UCN;
2494	stringstream ss;
2495	ss << UCN;
2496	string UCNstr = ss.str();
2497	string line;
2498	string strGcBasisLength;
2499	string strMonom, strCoeff, strCoeffNom, strCoeffDenom;
2500	int gcBasisLength=0;
2501	int intCoeff=1;
2502	int intCoeffNom=1; //better (number) but that's not compatible with stringstream;
2503	int intCoeffDenom=1;
2504	number nCoeff=nInit(1);
2505	number nCoeffNom=nInit(1);
2506	number nCoeffDenom=nInit(1);
2507	bool hasCoeffInQ = FALSE; //does the polynomial have rational coeff?
2508	bool hasNegCoeff = FALSE; //or a negative one?
2509	size_t found; //used for find_first_(not)_of
2510
2511	string prefix="/tmp/cone";
2512	string suffix=".sg";
2513	string filename;
2514	filename.append(prefix);
2515	filename.append(UCNstr);
2516	filename.append(suffix);
2517
2518	ifstream gcInputFile(filename.c_str(), ifstream::in);
2519
2520	ring saveRing=currRing;
2521	rChangeCurrRing(gc->baseRing);
2522	poly strPoly=pInit();
2523	poly resPoly=pInit(); //The poly to be read in
2524
2525	string::iterator EOL;
2526	int terms=1; //#Terms in the poly
2527
2528	while( !gcInputFile.eof() )
2529	{
2530	getline(gcInputFile,line);
2531	hasCoeffInQ = FALSE;
2532	hasNegCoeff = FALSE;
2533
2534	if(line=="GCBASISLENGTH")
2535	{
2536	getline(gcInputFile, line);
2537	strGcBasisLength = line;
2538	gcBasisLength=atoi(strGcBasisLength.c_str());
2539	}
2540	if(line=="GCBASIS")
2541	{
2542	for(int jj=0;jj<gcBasisLength;jj++)
2543	{
2544	getline(gcInputFile,line);
2545	//find out how many terms the poly consists of
2546	for(EOL=line.begin(); EOL!=line.end();++EOL)
2547	{
2548	string s;
2549	s=*EOL;
2550	if(s=="+" \|\| s=="-")
2551	terms++;
2552	}
2553	//magically convert strings into polynomials
2554	//polys.cc:p_Read
2555	//check until first occurance of + or -
2556	//data or c_str
2557	while(!line.empty())
2558	{
2559	hasNegCoeff = FALSE;
2560	found = line.find_first_of("+-"); //get the first monomial
2561	string tmp;
2562	tmp=line[found];
2563	// if(found!=0 && (tmp.compare("-")==0) )
2564	// hasNegCoeff = TRUE; //We have a coeff < 0
2565	if(found==0)
2566	{
2567	if(tmp.compare("-")==0)
2568	hasNegCoeff = TRUE;
2569	line.erase(0,1); //remove leading + or -
2570	found = line.find_first_of("+-"); //adjust found
2571	}
2572	strMonom = line.substr(0,found);
2573	line.erase(0,found);
2574	found = strMonom.find_first_of("/");
2575	if(found!=string::npos) //i.e. "/" exists in strMonom
2576	{
2577	hasCoeffInQ = TRUE;
2578	strCoeffNom=strMonom.substr(0,found);
2579	strCoeffDenom=strMonom.substr(found+1,strMonom.find_first_not_of("1234567890"));
2580	strMonom.erase(0,found);
2581	strMonom.erase(0,strMonom.find_first_not_of("1234567890/"));
2582	// ss << strCoeffNom;
2583	// ss >> intCoeffNom;
2584	// nCoeffNom=(snumber*)strCoeffNom.c_str();
2585	nRead(strCoeffNom.c_str(), &nCoeffNom);
2586	nRead(strCoeffDenom.c_str(), &nCoeffDenom);
2587	// nCoeffNom=nInit(intCoeffNom);
2588	// ss << strCoeffDenom;
2589	// ss >> intCoeffDenom;
2590	// nCoeffDenom=nInit(intCoeffDenom);
2591	// nCoeffDenom=(snumber*)strCoeffNom.c_str();
2592	//NOTE NOT SURE WHETHER THIS WILL WORK!
2593	//nCoeffNom=nInit(intCoeffNom);
2594	// nCoeffDenom=(number)strCoeffDenom.c_str();
2595	// nCoeffDenom=(number)strCoeffDenom.c_str();
2596	}
2597	else
2598	{
2599	found = strMonom.find_first_not_of("1234567890");
2600	strCoeff = strMonom.substr(0,found);
2601	if(!strCoeff.empty())
2602	{
2603	nRead(strCoeff.c_str(),&nCoeff);
2604	}
2605	else
2606	{
2607	// intCoeff = 1;
2608	nCoeff = nInit(1);
2609	}
2610
2611	}
2612	const char* monom = strMonom.c_str();
2613
2614	p_Read(monom,strPoly,currRing);
2615	switch (hasCoeffInQ)
2616	{
2617	case TRUE:
2618	if(hasNegCoeff)
2619	nCoeffNom=nNeg(nCoeffNom);
2620	// intCoeffNom *= -1;
2621	// pSetCoeff(strPoly, nDiv((number)intCoeffNom, (number)intCoeffDenom));
2622	pSetCoeff(strPoly, nDiv(nCoeffNom, nCoeffDenom));
2623	break;
2624	case FALSE:
2625	if(hasNegCoeff)
2626	nCoeff=nNeg(nCoeff);
2627	// intCoeff *= -1;
2628	if(!nIsOne(nCoeff))
2629	{
2630	// if(hasNegCoeff)
2631	// intCoeff *= -1;
2632	// pSetCoeff(strPoly,(number) intCoeff);
2633	pSetCoeff(strPoly, nCoeff );
2634	}
2635	break;
2636
2637	}
2638	//pSetCoeff(strPoly, (number) intCoeff);//Why is this set to zero instead of 1???
2639	if(pIsConstantComp(resPoly))
2640	resPoly=pCopy(strPoly);
2641	else
2642	resPoly=pAdd(resPoly,strPoly);
2643
2644
2645	}//while(!line.empty())
2646
2647	gcBasis->m[jj]=pCopy(resPoly);
2648	resPoly=NULL; //reset
2649	}
2650	break;
2651	}
2652	}//while(!gcInputFile.eof())
2653
2654	gcInputFile.close();
2655	rChangeCurrRing(saveRing);
2656	}
2657
2658
2659	int gcone::counter=0;
2660	int gfanHeuristic;
2661	ideal gfan(ideal inputIdeal, int h)
2662	{
2663	time_t tic,tac;
2664	time(&tic);
2665
2666	int numvar = pVariables;
2667	gfanHeuristic = h;
2668
2669	enum searchMethod {
2670	reverseSearch,
2671	noRevS
2672	};
2673
2674	searchMethod method;
2675	method = noRevS;
2676	// method = reverseSearch;
2677
2678	#ifdef gfan_DEBUG
2679	cout << "Now in subroutine gfan" << endl;
2680	#endif
2681	ring inputRing=currRing; // The ring the user entered
2682	ring rootRing; // The ring associated to the target ordering
2683	ideal res;
2684	facet *fRoot;
2685
2686	if (method==reverseSearch)
2687	{
2688
2689	/* Construct a new ring which will serve as our root*/
2690	rootRing=rCopy0(currRing);
2691	rootRing->order[0]=ringorder_lp;
2692
2693	rComplete(rootRing);
2694	rChangeCurrRing(rootRing);
2695
2696	/* Fetch the inputIdeal into our rootRing */
2697	map theMap=(map)idMaxIdeal(1); //evil hack!
2698	theMap->preimage=NULL; //neccessary?
2699	ideal rootIdeal;
2700	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
2701	#ifndef NDEBUG
2702	#ifdef gfan_DEBUG
2703	cout << "Root ideal is " << endl;
2704	idShow(rootIdeal);
2705	cout << "The root ring is " << endl;
2706	rWrite(rootRing);
2707	cout << endl;
2708	#endif
2709	#endif
2710
2711	//gcone *gcRoot = new gcone(); //Instantiate the sink
2712	gcone *gcRoot = new gcone(rootRing,rootIdeal);
2713	gcone *gcAct;
2714	gcAct = gcRoot;
2715	gcAct->numVars=pVariables;
2716	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
2717	#ifndef NDEBUG
2718	idShow(gcAct->gcBasis);
2719	#endif
2720	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
2721	//gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
2722	/*Now it is time to compute the search facets, respectively start the reverse search.
2723	But since we are in the root all facets should be search facets. IS THIS TRUE?
2724	NOTE: Check for flippability is not very sophisticated
2725	*/
2726	//gcAct->reverseSearch(gcAct);
2727	rChangeCurrRing(rootRing);
2728	res=gcRoot->gcBasis;
2729	}//if method==reverSearch
2730
2731	if(method==noRevS)
2732	{
2733	gcone *gcRoot = new gcone(currRing,inputIdeal);
2734	gcone *gcAct;
2735	gcAct = gcRoot;
2736	gcAct->numVars=pVariables;
2737	gcAct->getGB(inputIdeal);
2738	cout << "GB of input ideal is:" << endl;
2739	#ifndef NDEBUG
2740	idShow(gcAct->gcBasis);
2741	#endif
2742	if(gcAct->isMonomial(gcAct->gcBasis))
2743	{
2744	WerrorS("Monomial input - terminating");
2745	res=gcAct->gcBasis;
2746	//break;
2747	}
2748	cout << endl;
2749	gcAct->getConeNormals(gcAct->gcBasis);
2750	gcAct->noRevS(*gcAct);
2751	//res=gcAct->gcBasis;
2752	//Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
2753	res = inputIdeal;
2754	}
2755	dd_free_global_constants();
2756	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
2757	The return type will then be of type LIST_CMD
2758	Assume gfan has finished, thus we have enumerated all the cones
2759	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
2760	Groebner Basis and merge this somehow with LIST_CMD
2761	=> Count the cones!
2762	*/
2763	//rChangeCurrRing(rootRing);
2764	//res=gcAct->gcBasis;
2765	//res=gcRoot->gcBasis;
2766	time(&tac);
2767	cout << "Time: " << difftime(tac,tic) << "sec" << endl;
2768	return res;
2769	//return GBlist;
2770	}
2771	/*
2772	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
2773	*/
2774	#endif

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

Download in other formats: