Context Navigation

source: git/kernel/gfan.cc @ bb503c7

spielwiese

Last change on this file since bb503c7 was bb503c7, checked in by Martin Monerjan, 15 years ago
Cleanup idNorm in gcone::getGB git-svn-id: file:///usr/local/Singular/svn/trunk@12197 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.0 KB

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

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

Download in other formats: