Context Navigation

source: git/kernel/gfan.cc @ 7d60f7

spielwiese

Last change on this file since 7d60f7 was 7d60f7, checked in by Martin Monerjan, 14 years ago
kNF and timings git-svn-id: file:///usr/local/Singular/svn/trunk@12216 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.4 KB

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

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

Download in other formats: