Context Navigation

source: git/kernel/gfan.cc @ c5940e

spielwiese

Last change on this file since c5940e was c5940e, checked in by Martin Monerjan, 15 years ago
gcones doubly linked gcone::~gcone deletion of ideals in gcRoot-gcTmp list git-svn-id: file:///usr/local/Singular/svn/trunk@12149 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 67.8 KB

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

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

Download in other formats: