Context Navigation

source: git/kernel/gfan.cc @ 913422

spielwiese

Last change on this file since 913422 was 913422, checked in by Martin Monerjan, 14 years ago
gcone::isMonomial. Hack to avoid infinite loop when inputGB is monomial git-svn-id: file:///usr/local/Singular/svn/trunk@12201 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.3 KB

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

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

Download in other formats: