Context Navigation

source: git/kernel/gfan.cc @ d7ce08d

spielwiese

Last change on this file since d7ce08d was d7ce08d, checked in by Martin Monerjan, 15 years ago
Fixed normalize git-svn-id: file:///usr/local/Singular/svn/trunk@11849 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 48.3 KB

Line
1	/*
2	Compute the Groebner fan of an ideal
3	$Author: monerjan $
4	$Date: 2009-05-28 06:27:09 $
5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.56 2009-05-28 06:27:09 monerjan Exp $
6	$Id: gfan.cc,v 1.56 2009-05-28 06:27:09 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 <gmpxx.h>
28
29	/DO NOT REMOVE THIS/
30	#ifndef GMPRATIONAL
31	#define GMPRATIONAL
32	#endif
33
34	//Hacks for different working places
35	#define ITWM
36
37	#ifdef UNI
38	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
39	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
40	#endif
41
42	#ifdef HOME
43	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
44	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
45	#endif
46
47	#ifdef ITWM
48	#include "/u/slg/monerjan/cddlib/include/setoper.h"
49	#include "/u/slg/monerjan/cddlib/include/cdd.h"
50	#include "/u/slg/monerjan/cddlib/include/cddmp.h"
51	#endif
52
53	#ifndef gfan_DEBUG
54	#define gfan_DEBUG
55	#endif
56
57	//#include gcone.h
58	using namespace std;
59
60	/**
61	*\brief Class facet
62	* Implements the facet structure as a linked list
63	*
64	*/
65	class facet
66	{
67	private:
68	/** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
69	intvec *fNormal;
70
71	/** \brief An interior point of the facet*/
72	intvec *interiorPoint;
73
74	/** \brief Universal Cone Number
75	* The number of the cone the facet belongs to
76	*/
77	int UCN;
78
79	/** \brief The Groebner basis on the other side of a shared facet
80	*
81	* In order not to have to compute the flipped GB twice we store the basis we already get
82	* when identifying search facets. Thus in the next step of the reverse search we can
83	* just copy the old cone and update the facet and the gcBasis.
84	* facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
85	*/
86	ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip
87
88	public:
89	//bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i]
90	bool isIncoming; //Is the facet incoming or outgoing?
91	facet *next; //Pointer to next facet
92	//intvec codim2Normals =(intvec)omAlloc0(sizeof(intvec*)); //Integer matrix containing the (codim-2)-facets
93	struct c2N{
94	intvec normal;
95	intvec *next;
96	};
97
98	/** The default constructor. Do I need a constructor of type facet(intvec)? */
99	facet()
100	{
101	// Pointer to next facet. */
102	/* Defaults to NULL. This way there is no need to check explicitly */
103	this->next=NULL;
104	this->UCN=NULL;
105	//this->codim2Normals=NULL;
106	}
107
108	/** The default destructor */
109	~facet(){;}
110
111	/** Stores the facet normal \param intvec*/
112	void setFacetNormal(intvec *iv)
113	{
114	this->fNormal = ivCopy(iv);
115	//return;
116	}
117
118	/** Hopefully returns the facet normal */
119	intvec *getFacetNormal()
120	{
121	//this->fNormal->show(); cout << endl;
122	return this->fNormal;
123	}
124
125	/** Method to print the facet normal*/
126	void printNormal()
127	{
128	fNormal->show();
129	}
130
131	/** Store the flipped GB*/
132	void setFlipGB(ideal I)
133	{
134	this->flipGB=I;
135	}
136
137	/** Return the flipped GB*/
138	ideal getFlipGB()
139	{
140	return this->flipGB;
141	}
142
143	/** Print the flipped GB*/
144	void printFlipGB()
145	{
146	idShow(this->flipGB);
147	}
148
149	void setUCN(int n)
150	{
151	this->UCN=n;
152	}
153
154	int getUCN()
155	{
156	return this->UCN;
157	}
158
159	void setInteriorPoint(intvec *iv)
160	{
161	this->interiorPoint = ivCopy(iv);
162	}
163
164	intvec *getInteriorPoint()
165	{
166	return this->interiorPoint;
167	}
168
169	/*bool isFlippable(intvec &load)
170	{
171	bool res=TRUE;
172	int jj;
173	for (int jj = 0; jj<load.length(); jj++)
174	{
175	intvec *ivCanonical = new intvec(load.length());
176	(*ivCanonical)[jj]=1;
177	if (ivMult(&load,ivCanonical)<0)
178	{
179	res=FALSE;
180	break;
181	}
182	}
183	return res;
184
185	/*while (dotProduct(load,ivCanonical)>=0)
186	{
187	if (jj!=this->numVars)
188	{
189	intvec *ivCanonical = new intvec(this->numVars);
190	(*ivCanonical)[jj]=1;
191	res=TRUE;
192	jj += 1;
193	}
194	}
195	if (jj==this->numVars)
196	{
197	delete ivCanonical;
198	return FALSE;
199	}
200	else
201	{
202	delete ivCanonical;
203	return TRUE;
204	}*/
205	//}//bool isFlippable(facet &f)
206
207
208	friend class gcone; //Bad style
209	};
210
211	/**
212	*\brief Implements the cone structure
213	*
214	* A cone is represented by a linked list of facet normals
215	* @see facet
216	*/
217	/*class gcone
218	finally this should become s.th. like gconelib.{h,cc} to provide an API
219	*/
220	class gcone
221	{
222	private:
223	ring rootRing; //good to know this -> generic walk
224	ideal inputIdeal; //the original
225	ring baseRing; //the basering of the cone
226	/* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
227	intvec *ivIntPt; //an interior point of the cone
228	int UCN; //unique number of the cone
229
230	public:
231	/** \brief Pointer to the first facet */
232	facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals
233
234	/** # of variables in the ring */
235	int numVars; //#of variables in the ring
236
237	/** # of facets of the cone
238	* This value is set by gcone::getConeNormals
239	*/
240	int numFacets; //#of facets of the cone
241
242	/** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
243	ideal gcBasis; //GB of the cone, set by gcone::getGB();
244	gcone next; //Pointer to previous* cone in search tree
245
246	/** \brief Default constructor.
247	*
248	* Initialises this->next=NULL and this->facetPtr=NULL
249	*/
250	gcone()
251	{
252	this->next=NULL;
253	this->facetPtr=NULL;
254	this->baseRing=currRing;
255	this->UCN=1;
256	}
257
258	/** \brief Constructor with ring and ideal
259	*
260	* This constructor takes the root ring and the root ideal as parameters and stores
261	* them in the private members gcone::rootRing and gcone::inputIdeal
262	* Since knowledge of the root ring is only needed when using reverse search,
263	* this constructor is not needed when using the "second" method
264	*/
265	gcone(ring r, ideal I)
266	{
267	this->next=NULL;
268	this->facetPtr=NULL;
269	this->rootRing=r;
270	this->inputIdeal=I;
271	this->baseRing=currRing;
272	this->UCN=1;
273	}
274
275	/** \brief Copy constructor
276	*
277	* Copies one cone, sets this->gcBasis to the flipped GB and reverses the
278	* direction of the according facet normal
279	*/
280	//NOTE Prolly need to specify the facet to flip over
281	gcone(const gcone& gc)
282	{
283	this->next=NULL;
284	this->numVars=gc.numVars;
285	this->UCN=(gc.UCN)+1; //add 1 to the UCN of previous cone
286	facet *fAct= new facet();
287	this->facetPtr=fAct;
288
289	intvec *ivtmp = new intvec(this->numVars);
290	ivtmp = gc.facetPtr->getFacetNormal();
291
292	ideal gb;
293	gb=gc.facetPtr->getFlipGB();
294	this->gcBasis=gb; //this cone's GB is the flipped GB
295
296	/Reverse direction of the facet normal to make it an inner normal/
297	for (int ii=0; ii<this->numVars;ii++)
298	{
299	(ivtmp)[ii]=-(ivtmp)[ii];
300	}
301
302	fAct->setFacetNormal(ivtmp);
303	delete ivtmp;
304	delete fAct;
305	}
306
307	/** \brief Default destructor */
308	~gcone(){;} //destructor
309
310
311	/** \brief Set the interior point of a cone */
312	void setIntPoint(intvec *iv)
313	{
314	this->ivIntPt=ivCopy(iv);
315	}
316
317	/** \brief Return the interior point */
318	intvec *getIntPoint()
319	{
320	return this->ivIntPt;
321	}
322
323	void showIntPoint()
324	{
325	ivIntPt->show();
326	}
327
328	void showFacets()
329	{
330	facet *f = new facet();
331	f = this->facetPtr;
332	intvec *iv = new intvec(this->numVars);
333
334	while (f->next!=NULL)
335	{
336	iv = f->getFacetNormal();
337	iv->show();
338	f=f->next;
339	}
340	//delete iv;
341	//delete f;
342	}
343
344	/** \brief Set gcone::numFacets */
345	void setNumFacets()
346	{
347	}
348
349	/** \brief Get gcone::numFacets */
350	int getNumFacets()
351	{
352	return this->numFacets;
353	}
354
355	int getUCN()
356	{
357	return this->UCN;
358	}
359
360	/** \brief Compute the normals of the cone
361	*
362	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
363	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
364	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
365	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
366	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
367	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
368	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
369	*
370	* Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
371	* an interior point of the cone.
372	*/
373	void getConeNormals(ideal const &I, bool compIntPoint=FALSE)
374	{
375	#ifdef gfan_DEBUG
376	std::cout << "* Computing Inequalities... *" << std::endl;
377	#endif
378	//All variables go here - except ineq matrix and *v, see below
379	int lengthGB=IDELEMS(I); // # of polys in the groebner basis
380	int pCompCount; // # of terms in a poly
381	poly aktpoly;
382	int numvar = pVariables; // # of variables in a polynomial (or ring?)
383	int leadexp[numvar]; // dirty hack of exp.vects
384	int aktexp[numvar];
385	int cols,rows; // will contain the dimensions of the ineq matrix - deprecated by
386	dd_rowrange ddrows;
387	dd_colrange ddcols;
388	dd_rowset ddredrows; // # of redundant rows in ddineq
389	dd_rowset ddlinset; // the opposite
390	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
391	dd_NumberType ddnumb=dd_Integer; //Number type
392	dd_ErrorType dderr=dd_NoError; //
393	// End of var declaration
394	#ifdef gfan_DEBUG
395	cout << "The Groebner basis has " << lengthGB << " elements" << endl;
396	cout << "The current ring has " << numvar << " variables" << endl;
397	#endif
398	cols = numvar;
399
400	//Compute the # inequalities i.e. rows of the matrix
401	rows=0; //Initialization
402	for (int ii=0;ii<IDELEMS(I);ii++)
403	{
404	aktpoly=(poly)I->m[ii];
405	rows=rows+pLength(aktpoly)-1;
406	}
407	#ifdef gfan_DEBUG
408	cout << "rows=" << rows << endl;
409	cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
410	#endif
411	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
412	dd_set_global_constants();
413	ddrows=rows;
414	ddcols=cols;
415	dd_MatrixPtr ddineq; //Matrix to store the inequalities
416	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
417
418	// We loop through each g\in GB and compute the resulting inequalities
419	for (int i=0; i<IDELEMS(I); i++)
420	{
421	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
422	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
423	#ifdef gfan_DEBUG
424	cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
425	#endif
426
427	int v=(int )omAlloc((numvar+1)*sizeof(int));
428	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
429
430	//Store leadexp for aktpoly
431	for (int kk=0;kk<numvar;kk++)
432	{
433	leadexp[kk]=v[kk+1];
434	//Since we need to know the difference of leadexp with the other expvects we do nothing here
435	//but compute the diff below
436	}
437
438
439	while (pNext(aktpoly)!=NULL) //move to next term until NULL
440	{
441	aktpoly=pNext(aktpoly);
442	pSetm(aktpoly); //doesn't seem to help anything
443	pGetExpV(aktpoly,v);
444	for (int kk=0;kk<numvar;kk++)
445	{
446	aktexp[kk]=v[kk+1];
447	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
448	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
449	}
450	aktmatrixrow=aktmatrixrow+1;
451	}//while
452
453	} //for
454
455	//Maybe add another row to contain the constraints of the standard simplex?
456
457	#ifdef gfan_DEBUG
458	cout << "The inequality matrix is" << endl;
459	dd_WriteMatrix(stdout, ddineq);
460	#endif
461
462	// The inequalities are now stored in ddineq
463	// Next we check for superflous rows
464	ddredrows = dd_RedundantRows(ddineq, &dderr);
465	if (dderr!=dd_NoError) // did an error occur?
466	{
467	dd_WriteErrorMessages(stderr,dderr); //if so tell us
468	} else
469	{
470	cout << "Redundant rows: ";
471	set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
472	}//if dd_Error
473
474	//Remove reduntant rows here!
475	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
476	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
477	ddcols = ddineq->colsize;
478	#ifdef gfan_DEBUG
479	cout << "Having removed redundancies, the normals now read:" << endl;
480	dd_WriteMatrix(stdout,ddineq);
481	cout << "Rows = " << ddrows << endl;
482	cout << "Cols = " << ddcols << endl;
483	#endif
484
485	/Write the normals into class facet/
486	#ifdef gfan_DEBUG
487	cout << "Creating list of normals" << endl;
488	#endif
489	/The pointer fRoot should be the return value of this function*/
490	facet *fRoot = new facet(); //instantiate new facet
491	this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
492	facet *fAct; //instantiate pointer to active facet
493	fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress!
494	//std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
495	for (int kk = 0; kk<ddrows; kk++)
496	{
497	intvec *load = new intvec(this->numVars); //intvec to store a single facet normal that will then be stored via setFacetNormal
498	for (int jj = 1; jj <ddcols; jj++)
499	{
500	double foo;
501	foo = mpq_get_d(ddineq->matrix[kk][jj]);
502	#ifdef gfan_DEBUG
503	std::cout << "fAct is " << foo << " at " << fAct << std::endl;
504	#endif
505	(*load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
506	}//for (int jj = 1; jj <ddcols; jj++)
507
508	/Quick'n'dirty hack for flippability/
509	bool isFlippable=FALSE;
510	for (int jj = 0; jj<load->length(); jj++)
511	{
512	intvec *ivCanonical = new intvec(load->length());
513	(*ivCanonical)[jj]=1;
514	cout << "dotProd=" << dotProduct(load,ivCanonical) << endl;
515	if (dotProduct(load,ivCanonical)<0)
516	//if (ivMult(load,ivCanonical)<0)
517	{
518	isFlippable=TRUE;
519	break; //URGHS
520	}
521	}
522	if (isFlippable==FALSE)
523	{
524	cout << "Ignoring facet";
525	load->show();
526	//fAct->next=NULL;
527	}
528	else
529	{ /Now load should be full and we can call setFacetNormal/
530	fAct->setFacetNormal(load);
531	fAct->next = new facet();
532	fAct=fAct->next; //this should definitely not be called in the above while-loop :D
533	}//if (isFlippable==FALSE)
534	//delete load;
535	}//for (int kk = 0; kk<ddrows; kk++)
536
537	/*
538	Now we should have a linked list containing the facet normals of those facets that are
539	-irredundant
540	-flipable
541	Adressing is done via *facetPtr
542	*/
543
544	if (compIntPoint==TRUE)
545	{
546	intvec *iv = new intvec(this->numVars);
547	interiorPoint(ddineq, *iv); //NOTE ddineq contains non-flippable facets
548	this->setIntPoint(iv); //stores the interior point in gcone::ivIntPt
549	//delete iv;
550	}
551
552	//Compute the number of facets
553	this->numFacets=ddineq->rowsize;
554
555	//Clean up but don't delete the return value! (Whatever it will turn out to be)
556	//dd_FreeMatrix(ddineq);
557	//set_free(ddredrows);
558	//free(ddnewpos);
559	//set_free(ddlinset);
560	//NOTE Commented out. Solved the bug that after facet2Matrix there were facets lost
561	//THIS SUCKS
562	//dd_free_global_constants();
563
564	}//method getConeNormals(ideal I)
565
566
567	/** \brief Compute the Groebner Basis on the other side of a shared facet
568	*
569	* Implements algorithm 4.3.2 from Anders' thesis.
570	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
571	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
572	* is parallel to \f$ leadexp(g) \f$
573	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
574	* Other possibilities include computing the rank of the matrix consisting of the vectors in question and
575	* computing an interior point of the facet and taking all terms having the same weight with respect
576	* to this interior point.
577	*\param ideal, facet
578	* Input: a marked,reduced Groebner basis and a facet
579	*/
580	void flip(ideal gb, facet *f) //Compute "the other side"
581	{
582	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
583	fNormal = f->getFacetNormal(); //read this->fNormal;
584	#ifdef gfan_DEBUG
585	std::cout << "===" << std::endl;
586	std::cout << "running gcone::flip" << std::endl;
587	// std::cout << "fNormal=";
588	// fNormal->show();
589	// std::cout << std::endl;
590	#endif
591	/1st step: Compute the initial ideal/
592	poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form
593	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
594	poly aktpoly;
595	intvec *check = new intvec(this->numVars); //array to store the difference of LE and v
596
597	for (int ii=0;ii<IDELEMS(gb);ii++)
598	{
599	aktpoly = (poly)gb->m[ii];
600	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
601	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
602	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
603	initialFormElement[ii]=pHead(aktpoly);
604
605	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
606	{
607	aktpoly=pNext(aktpoly); //next term
608	pSetm(aktpoly);
609	pGetExpV(aktpoly,v);
610	/* Convert (int)v into (intvec)check */
611	for (int jj=0;jj<this->numVars;jj++)
612	{
613	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
614	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
615	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
616	}
617	#ifdef gfan_DEBUG
618	// cout << "check=";
619	// check->show();
620	// cout << endl;
621	#endif
622	//TODO why not check, fNormal????
623	if (isParallel(check,fNormal)) //pass *check when
624	{
625	// cout << "Parallel vector found, adding to initialFormElement" << endl;
626	initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
627	}
628	}//while
629	#ifdef gfan_DEBUG
630	cout << "Initial Form=";
631	pWrite(initialFormElement[ii]);
632	cout << "---" << endl;
633	#endif
634	/Now initialFormElement must be added to (ideal)initialForm /
635	initialForm->m[ii]=initialFormElement[ii];
636	}//for
637	#ifdef gfan_DEBUG
638	cout << "Initial ideal is: " << endl;
639	idShow(initialForm);
640	//f->printFlipGB();
641	cout << "===" << endl;
642	#endif
643	//delete check;
644
645	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
646	/*Substep 2.1
647	compute $G_{-\alpha}(in_v(I))
648	see journal p. 66
649	*/
650	ring srcRing=currRing;
651
652	//intvec *negfNormal = new intvec(this->numVars);
653	//negfNormal=ivNeg(fNormal);
654	ring tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
655	rChangeCurrRing(tmpRing);
656
657	//rWrite(currRing); cout << endl;
658
659	ideal ina;
660	ina=idrCopyR(initialForm,srcRing);
661	#ifdef gfan_DEBUG
662	cout << "ina=";
663	idShow(ina); cout << endl;
664	#endif
665	ideal H;
666	H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
667	idSkipZeroes(H);
668	#ifdef gfan_DEBUG
669	// cout << "H="; idShow(H); cout << endl;
670	#endif
671	/*Substep 2.2
672	do the lifting and mark according to H
673	*/
674	rChangeCurrRing(srcRing);
675	ideal srcRing_H;
676	ideal srcRing_HH;
677	srcRing_H=idrCopyR(H,tmpRing);
678	#ifdef gfan_DEBUG
679	// cout << "srcRing_H = ";
680	// idShow(srcRing_H); cout << endl;
681	#endif
682	srcRing_HH=ffG(srcRing_H,this->gcBasis);
683	#ifdef gfan_DEBUG
684	// cout << "srcRing_HH = ";
685	// idShow(srcRing_HH); cout << endl;
686	#endif
687	/*Substep 2.2.1
688	Mark according to G_-\alpha
689	Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
690	we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
691	represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
692	Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
693	compute the difference accordingly
694	*/
695	dd_set_global_constants();
696	bool markingsAreCorrect=FALSE;
697	dd_MatrixPtr intPointMatrix;
698	int iPMatrixRows=0;
699	dd_rowrange aktrow=0;
700	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
701	{
702	poly aktpoly=(poly)srcRing_HH->m[ii];
703	iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
704	}
705	/* additionally one row for the standard-simplex and another for a row that becomes 0 during
706	construction of the differences
707	*/
708	intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
709	intPointMatrix->numbtype=dd_Integer; //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
710
711	for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
712	{
713	markingsAreCorrect=FALSE; //crucial to initialise here
714	poly aktpoly=srcRing_HH->m[ii];
715	/Comparison of leading monomials is done via exponent vectors/
716	for (int jj=0;jj<IDELEMS(H);jj++)
717	{
718	int src_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
719	int dst_ExpV = (int )omAlloc((this->numVars+1)*sizeof(int));
720	pGetExpV(aktpoly,src_ExpV);
721	rChangeCurrRing(tmpRing); //this ring change is crucial!
722	pGetExpV(pCopy(H->m[ii]),dst_ExpV);
723	rChangeCurrRing(srcRing);
724	bool expVAreEqual=TRUE;
725	for (int kk=1;kk<=this->numVars;kk++)
726	{
727	#ifdef gfan_DEBUG
728	//cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
729	#endif
730	if (src_ExpV[kk]!=dst_ExpV[kk])
731	{
732	expVAreEqual=FALSE;
733	}
734	}
735	//if (src_ExpV == dst_ExpV)
736	if (expVAreEqual==TRUE)
737	{
738	markingsAreCorrect=TRUE; //everything is fine
739	#ifdef gfan_DEBUG
740	// cout << "correct markings" << endl;
741	#endif
742	}//if (pHead(aktpoly)==pHead(H->m[jj])
743	delete src_ExpV;
744	delete dst_ExpV;
745	}//for (int jj=0;jj<IDELEMS(H);jj++)
746
747	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
748	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
749	if (markingsAreCorrect==TRUE)
750	{
751	pGetExpV(aktpoly,leadExpV);
752	}
753	else
754	{
755	rChangeCurrRing(tmpRing);
756	pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
757	rChangeCurrRing(srcRing);
758	}
759	/compute differences of the expvects/
760	while (pNext(aktpoly)!=NULL)
761	{
762	/*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
763	is not omitted when computing the differences*/
764	if(markingsAreCorrect==TRUE)
765	{
766	aktpoly=pNext(aktpoly);
767	pGetExpV(aktpoly,v);
768	}
769	else
770	{
771	pGetExpV(pHead(aktpoly),v);
772	markingsAreCorrect=TRUE;
773	}
774
775	for (int jj=0;jj<this->numVars;jj++)
776	{
777	/Store into ddMatrix/
778	dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
779	}
780	aktrow +=1;
781	}
782	delete v;
783	delete leadExpV;
784	}//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
785	/Now we add the constraint for the standard simplex/
786	/NOTE:Might actually work without the standard simplex/
787	dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
788	for (int jj=1;jj<=this->numVars;jj++)
789	{
790	dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
791	}
792	dd_WriteMatrix(stdout,intPointMatrix);
793	intvec *iv_weight = new intvec(this->numVars);
794	interiorPoint(intPointMatrix, *iv_weight); //iv_weight now contains the interior point
795	dd_FreeMatrix(intPointMatrix);
796	dd_free_global_constants();
797
798	/*Step 3
799	turn the minimal basis into a reduced one
800	*/
801	int i,j;
802	ring dstRing=rCopy0(srcRing);
803	i=rBlocks(srcRing);
804
805	dstRing->order=(int )omAlloc((i+1)sizeof(int));
806	for(j=i;j>0;j--)
807	{
808	dstRing->order[j]=srcRing->order[j-1];
809	dstRing->block0[j]=srcRing->block0[j-1];
810	dstRing->block1[j]=srcRing->block1[j-1];
811	if (srcRing->wvhdl[j-1] != NULL)
812	{
813	dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
814	}
815	}
816	dstRing->order[0]=ringorder_a;
817	dstRing->order[1]=ringorder_dp;
818	dstRing->order[2]=ringorder_C;
819	dstRing->wvhdl[0] =( int )omAlloc((iv_weight->length())sizeof(int));
820
821	for (int ii=0;ii<this->numVars;ii++)
822	{
823	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
824	}
825	rComplete(dstRing);
826
827	// NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
828	// Thus:
829	//ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
830	//cout << "PLING" << endl;
831	/*ring dstRing=rCopy0(srcRing);
832	for (int ii=0;ii<this->numVars;ii++)
833	{
834	dstRing->wvhdl[0][ii]=(*iv_weight)[ii];
835	}*/
836	rChangeCurrRing(dstRing);
837	#ifdef gfan_DEBUG
838	rWrite(dstRing); cout << endl;
839	#endif
840	ideal dstRing_I;
841	dstRing_I=idrCopyR(srcRing_HH,srcRing);
842	//validOpts<1>=TRUE;
843	#ifdef gfan_DEBUG
844	//idShow(dstRing_I);
845	#endif
846	BITSET save=test;
847	test\|=Sy_bit(OPT_REDSB);
848	test\|=Sy_bit(6); //OPT_DEBUG
849	dstRing_I=kStd(idrCopyR(this->inputIdeal,this->rootRing),NULL,testHomog,NULL);
850	kInterRed(dstRing_I);
851	idSkipZeroes(dstRing_I);
852	test=save;
853	/End of step 3 - reduction/
854
855	f->setFlipGB(dstRing_I);//store the flipped GB
856	#ifdef gfan_DEBUG
857	cout << "Flipped GB is: " << endl;
858	f->printFlipGB();
859	#endif
860	}//void flip(ideal gb, facet *f)
861
862	/** \brief Compute the remainder of a polynomial by a given ideal
863	*
864	* Compute \f$ f^{\mathcal{G}} \f$
865	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
866	* However, since we are only interested in the remainder, there is no need to
867	* compute the factors \f$ a_i \f$
868	*/
869	//NOTE: Should be replaced by kNF or kNF2
870	poly restOfDiv(poly const &f, ideal const &I)
871	{
872	cout << "Entering restOfDiv" << endl;
873	poly p=f;
874	//pWrite(p);
875	//poly r=kCreateZeroPoly(,currRing,currRing); //The 0-polynomial, hopefully
876	poly r=NULL; //The zero polynomial
877	int ii;
878	bool divOccured;
879
880	while (p!=NULL)
881	{
882	ii=1;
883	divOccured=FALSE;
884
885	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
886	{
887	if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
888	{
889	poly step1,step2,step3;
890	//cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
891	step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
892	//cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
893	step2 = ppMult_qq(step1, I->m[ii-1]);
894	step3 = pSub(pCopy(p), step2);
895	//p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
896	//pSetm(p);
897	pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
898	p=step3;
899	//pWrite(p);
900	divOccured=TRUE;
901	}
902	else
903	{
904	ii += 1;
905	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
906	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
907	if (divOccured==FALSE)
908	{
909	//cout << "TICK 5" << endl;
910	r=pAdd(pCopy(r),pHead(p));
911	pSetm(r);
912	pSort(r);
913	//cout << "r="; pWrite(r); cout << endl;
914
915	if (pLength(p)!=1)
916	{
917	p=pSub(pCopy(p),pHead(p)); //Here it may occur that p=0 instead of p=NULL
918	}
919	else
920	{
921	p=NULL; //Hack to correct this situation
922	}
923	//cout << "p="; pWrite(p);
924	}//if (divOccured==FALSE)
925	}//while (p!=0)
926	return r;
927	}//poly restOfDiv(poly const &f, ideal const &I)
928
929	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
930	*/
931	//NOTE: use kNF or kNF2 instead of restOfDivision
932	ideal ffG(ideal const &H, ideal const &G)
933	{
934	cout << "Entering ffG" << endl;
935	int size=IDELEMS(H);
936	ideal res=idInit(size,1);
937	poly temp1, temp2, temp3; //polys to temporarily store values for pSub
938	for (int ii=0;ii<size;ii++)
939	{
940	res->m[ii]=restOfDiv(H->m[ii],G);
941	//res->m[ii]=kNF(H->m[ii],G);
942	temp1=H->m[ii];
943	temp2=res->m[ii];
944	temp3=pSub(temp1, temp2);
945	res->m[ii]=temp3;
946	//res->m[ii]=pSub(temp1,temp2); //buggy
947	//pSort(res->m[ii]);
948	//pSetm(res->m[ii]);
949	//cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);
950	}
951	return res;
952	}
953
954	/** \brief Compute a Groebner Basis
955	*
956	* Computes the Groebner basis and stores the result in gcone::gcBasis
957	*\param ideal
958	*\return void
959	*/
960	void getGB(ideal const &inputIdeal)
961	{
962	ideal gb;
963	gb=kStd(inputIdeal,NULL,testHomog,NULL);
964	idSkipZeroes(gb);
965	this->gcBasis=gb; //write the GB into gcBasis
966	}//void getGB
967
968	/** \brief The Generic Groebner Walk due to FJLT
969	* Needed for computing the search facet
970	*/
971	ideal GenGrbWlk(ideal, ideal)
972	{
973	}//GGW
974
975	/** \brief Compute the negative of a given intvec
976	*/
977	intvec *ivNeg(intvec const &iv)
978	{
979	intvec *res = new intvec(this->numVars);
980	for(int ii=0;ii<this->numVars;ii++)
981	{
982	res[ii] = -(int)iv[ii];
983	}
984	return res;
985	}
986
987
988	/** \brief Compute the dot product of two intvecs
989	*
990	*/
991	int dotProduct(intvec const &iva, intvec const &ivb)
992	{
993	//intvec iva=a;
994	//intvec ivb=b;
995	int res=0;
996	for (int i=0;i<this->numVars;i++)
997	{
998	res = res+(iva[i]*ivb[i]);
999	}
1000	return res;
1001	}//int dotProduct
1002
1003	/** \brief Check whether two intvecs are parallel
1004	*
1005	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
1006	*/
1007	bool isParallel(intvec const &a, intvec const &b)
1008	{
1009	int lhs,rhs;
1010	lhs=dotProduct(a,b)*dotProduct(a,b);
1011	rhs=dotProduct(a,a)*dotProduct(b,b);
1012	//cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
1013	if (lhs==rhs)
1014	{
1015	return TRUE;
1016	}
1017	else
1018	{
1019	return FALSE;
1020	}
1021	}//bool isParallel
1022
1023	/** \brief Compute an interior point of a given cone
1024	* Result will be written into intvec iv
1025	*/
1026	void interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
1027	{
1028	dd_LPPtr lp,lpInt;
1029	dd_ErrorType err=dd_NoError;
1030	dd_LPSolverType solver=dd_DualSimplex;
1031	dd_LPSolutionPtr lpSol=NULL;
1032	dd_rowset ddlinset,ddredrows; //needed for dd_FindRelativeInterior
1033	dd_rowindex ddnewpos;
1034	dd_NumberType numb;
1035	//M->representation=dd_Inequality;
1036	//M->objective-dd_LPMin; //Not sure whether this is needed
1037
1038	//NOTE: Make this n-dimensional!
1039	//dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1040
1041	//dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1042	//if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1043	//cout << "Tick 2" << endl;
1044	//dd_WriteMatrix(stdout,M);
1045
1046	lp=dd_Matrix2LP(M, &err);
1047	if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}
1048	if (lp==NULL){cout << "LP is NULL" << endl;}
1049	#ifdef gfan_DEBUG
1050	// dd_WriteLP(stdout,lp);
1051	#endif
1052
1053	lpInt=dd_MakeLPforInteriorFinding(lp);
1054	if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1055	#ifdef gfan_DEBUG
1056	// dd_WriteLP(stdout,lpInt);
1057	#endif
1058
1059	dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1060	if (err!=dd_NoError)
1061	{
1062	cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1063	dd_WriteErrorMessages(stdout, err);
1064	}
1065
1066	//dd_LPSolve(lpInt,solver,&err); //This will not result in a point from the relative interior
1067	if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1068	//cout << "Tick 5" << endl;
1069
1070	//lpSol=dd_CopyLPSolution(lpInt);
1071	if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}
1072	//cout << "Tick 6" << endl;
1073	#ifdef gfan_DEBUG
1074	cout << "Interior point: ";
1075	#endif
1076	for (int ii=1; ii<(lpSol->d)-1;ii++)
1077	{
1078	#ifdef gfan_DEBUG
1079	dd_WriteNumber(stdout,lpSol->sol[ii]);
1080	#endif
1081	/* NOTE This works only as long as gmp returns fractions with the same denominator*/
1082	(iv)[ii-1]=(int)mpz_get_d(mpq_numref(lpSol->sol[ii])); //looks evil, but does the trick
1083	}
1084	dd_FreeLPSolution(lpSol);
1085	dd_FreeLPData(lpInt);
1086	dd_FreeLPData(lp);
1087	set_free(ddlinset);
1088	set_free(ddredrows);
1089
1090	}//void interiorPoint(dd_MatrixPtr const &M)
1091
1092	/** \brief Copy a ring and add a weighted ordering in first place
1093	* Kudos to walkSupport.cc
1094	*/
1095	ring rCopyAndAddWeight(ring const &r, intvec const *ivw)
1096	{
1097	ring res=(ring)omAllocBin(ip_sring_bin);
1098	memcpy4(res,r,sizeof(ip_sring));
1099	res->VarOffset = NULL;
1100	res->ref=0;
1101
1102	if (r->algring!=NULL)
1103	r->algring->ref++;
1104	if (r->parameter!=NULL)
1105	{
1106	res->minpoly=nCopy(r->minpoly);
1107	int l=rPar(r);
1108	res->parameter=(char *)omAlloc(lsizeof(char_ptr));
1109	int i;
1110	for(i=0;i<rPar(r);i++)
1111	{
1112	res->parameter[i]=omStrDup(r->parameter[i]);
1113	}
1114	}
1115
1116	int i=rBlocks(r);
1117	int jj;
1118
1119	res->order =(int )omAlloc((i+1)sizeof(int));
1120	res->block0=(int )omAlloc((i+1)sizeof(int));
1121	res->block1=(int )omAlloc((i+1)sizeof(int));
1122	res->wvhdl =(int *)omAlloc((i+1)sizeof(int**));
1123	for(jj=0;jj<i;jj++)
1124	{
1125	if (r->wvhdl[jj] != NULL)
1126	{
1127	res->wvhdl[jj] = (int*) omMemDup(r->wvhdl[jj-1]);
1128	}
1129	else
1130	{
1131	res->wvhdl[jj+1]=NULL;
1132	}
1133	}
1134
1135	for (jj=0;jj<i;jj++)
1136	{
1137	res->order[jj+1]=r->order[jj];
1138	res->block0[jj+1]=r->block0[jj];
1139	res->block1[jj+1]=r->block1[jj];
1140	}
1141
1142	res->order[0]=ringorder_a;
1143	res->order[1]=ringorder_dp; //basically useless, since that should never be used
1144	int length=ivw->length();
1145	int A=(int )omAlloc(length*sizeof(int));
1146	for (jj=0;jj<length;jj++)
1147	{
1148	A[jj]=(*ivw)[jj];
1149	}
1150	res->wvhdl[0]=(int *)A;
1151	res->block0[0]=1;
1152	res->block1[0]=length;
1153
1154	res->names = (char *)omAlloc0(rVar(r) sizeof(char_ptr));
1155	for (i=rVar(res)-1;i>=0; i--)
1156	{
1157	res->names[i] = omStrDup(r->names[i]);
1158	}
1159	rComplete(res);
1160	return res;
1161	}//rCopyAndAdd
1162
1163	ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1164	{
1165	ring res=rCopy0(currRing);
1166	rComplete(res);
1167	rSetWeightVec(res,(int64*)ivw);
1168	//rChangeCurrRing(rnew);
1169	return res;
1170	}
1171
1172	/** \brief Checks whether a given facet is a search facet
1173	* Determines whether a given facet of a cone is the search facet of a neighbouring cone
1174	* This is done in the following way:
1175	* We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1176	* that is first crossed during the generic walk.
1177	* We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1178	* If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1179	*/
1180	bool isSearchFacet(gcone &gcTmp, facet *testfacet)
1181	{
1182	ring actRing=currRing;
1183	facet facetPtr=(facet)gcTmp.facetPtr;
1184	facet fMin=new facet(facetPtr); //Pointer to the "minimal" facet
1185	//facet *fMin = new facet(tmpcone.facetPtr);
1186	//fMin=tmpcone.facetPtr; //Initialise to first facet of tmpcone
1187	facet *fAct; //Ptr to alpha_i
1188	facet *fCmp; //Ptr to alpha_j
1189	fAct = fMin;
1190	fCmp = fMin->next;
1191
1192	rChangeCurrRing(this->rootRing); //because we compare the monomials in the rootring
1193	poly p=pInit();
1194	poly q=pInit();
1195	intvec *alpha_i = new intvec(this->numVars);
1196	intvec *alpha_j = new intvec(this->numVars);
1197	intvec *sigma = new intvec(this->numVars);
1198	sigma=gcTmp.getIntPoint();
1199
1200	int u=(int )omAlloc((this->numVars+1)*sizeof(int));
1201	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
1202	u[0]=0; v[0]=0;
1203	int weight1,weight2;
1204	while(fAct->next->next!=NULL) //NOTE this is ugly. Can it be done without fCmp?
1205	{
1206	/* Get alpha_i and alpha_{i+1} */
1207	alpha_i=fAct->getFacetNormal();
1208	alpha_j=fCmp->getFacetNormal();
1209	#ifdef gfan_DEBUG
1210	alpha_i->show();
1211	alpha_j->show();
1212	#endif
1213	/Compute the dot product of sigma and alpha_{i,j}/
1214	weight1=dotProduct(sigma,alpha_i);
1215	weight2=dotProduct(sigma,alpha_j);
1216	#ifdef gfan_DEBUG
1217	cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1218	#endif
1219	/Adjust alpha_i and alpha_i+1 accordingly/
1220	for(int ii=1;ii<=this->numVars;ii++)
1221	{
1222	u[ii]=weight1(alpha_i)[ii-1];
1223	v[ii]=weight2(alpha_j)[ii-1];
1224	}
1225
1226	/Now p_weight and q_weight need to be compared as exponent vectors/
1227	pSetCoeff0(p,nInit(1));
1228	pSetCoeff0(q,nInit(1));
1229	pSetExpV(p,u);
1230	pSetm(p);
1231	pSetExpV(q,v);
1232	pSetm(q);
1233	#ifdef gfan_DEBUG
1234	pWrite(p);pWrite(q);
1235	#endif
1236	/*We want to check whether x^p < x^q
1237	=> want to check for return value 1 */
1238	if (pLmCmp(p,q)==1) //i.e. x^q is smaller
1239	{
1240	fMin=fCmp;
1241	fAct=fMin;
1242	fCmp=fCmp->next;
1243	}
1244	else
1245	{
1246	//fAct=fAct->next;
1247	if(fCmp->next!=NULL)
1248	{
1249	fCmp=fCmp->next;
1250	}
1251	else
1252	{
1253	fAct=fAct->next;
1254	}
1255	}
1256	//fAct=fAct->next;
1257	}//while(fAct.facetPtr->next!=NULL)
1258	delete alpha_i,alpha_j,sigma;
1259
1260	/If testfacet was minimal then fMin should still point there /
1261
1262	//if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))
1263	#ifdef gfan_DEBUG
1264	cout << "Checking for parallelity" << endl <<" fMin is";
1265	fMin->printNormal();
1266	cout << "testfacet is ";
1267	testfacet->printNormal();
1268	cout << endl;
1269	#endif
1270	if (fMin==gcTmp.facetPtr)
1271	//if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
1272	//if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
1273	{
1274	cout << "Parallel" << endl;
1275	rChangeCurrRing(actRing);
1276	//delete alpha_min, test;
1277	return TRUE;
1278	}
1279	else
1280	{
1281	cout << "Not parallel" << endl;
1282	rChangeCurrRing(actRing);
1283	//delete alpha_min, test;
1284	return FALSE;
1285	}
1286	}//bool isSearchFacet
1287
1288	/** \brief Check for equality of two intvecs
1289	*/
1290	bool areEqual(intvec const &a, intvec const &b)
1291	{
1292	bool res=TRUE;
1293	for(int ii=0;ii<this->numVars;ii++)
1294	{
1295	if(a[ii]!=b[ii])
1296	{
1297	res=FALSE;
1298	break;
1299	}
1300	}
1301	return res;
1302	}
1303
1304	/** \brief The reverse search algorithm
1305	*/
1306	void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1307	{
1308	facet *fAct=new facet();
1309	fAct = gcAct->facetPtr;
1310
1311	while(fAct->next!=NULL) //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
1312	{
1313	cout << "==========================================================================================="<< endl;
1314	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1315	gcone gcTmp = new gcone(gcAct);
1316	//idShow(gcTmp->gcBasis);
1317	gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1318	#ifdef gfan_DEBUG
1319	facet *f = new facet();
1320	f=gcTmp->facetPtr;
1321	while(f->next!=NULL)
1322	{
1323	f->printNormal();
1324	f=f->next;
1325	}
1326	#endif
1327	gcTmp->showIntPoint();
1328	/recursive part goes gere/
1329	if (isSearchFacet(gcTmp,(facet)gcAct->facetPtr))
1330	{
1331	gcAct->next=gcTmp;
1332	cout << "PING"<< endl;
1333	reverseSearch(gcTmp);
1334	}
1335	else
1336	{
1337	delete gcTmp;
1338	/NOTE remove fAct from linked list. It's no longer needed/
1339	}
1340	/recursion ends/
1341	fAct = fAct->next;
1342	}//while(fAct->next!=NULL)
1343	}//reverseSearch
1344
1345	/** \brief The new method of Markwig and Jensen
1346	* Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
1347	* In order to represent a facet uniquely compute also the \f$(codim-2)\f$-facets and norm by the gcd of the components.
1348	* Keep a list of facets as a linked list containing an intvec and an integer matrix.
1349	* Since a \f$(codim-1)\f$-facet belongs to exactly two full dimensional cones, we remove a facet from the list as
1350	* soon as we compute this facet again. Comparison of facets is done by...
1351	*/
1352	void noRevS(gcone &gc, bool usingIntPoint=FALSE)
1353	{
1354	facet *fListPtr = new facet(); //The list containing ALL facets we come across
1355	facet *fAct = new facet();
1356	fAct = gc.facetPtr;
1357
1358	#ifdef gfan_DEBUG
1359	cout << "NoRevs" << endl;
1360	cout << "Facets are:" << endl;
1361	gc.showFacets();
1362	#endif
1363
1364	dd_set_global_constants();
1365	dd_rowrange ddrows;
1366	dd_colrange ddcols;
1367	ddrows=2; //Each facet is described by its normal
1368	ddcols=gc.numVars+1; // plus one for the standard simplex
1369	if(usingIntPoint){
1370	while(fAct->next!=NULL)
1371	{
1372	/Compute an interior point for each facet/
1373	dd_MatrixPtr ddineq;
1374	ddineq=dd_CreateMatrix(ddrows,ddcols);
1375	intvec *comp = new intvec(this->numVars);
1376	comp=fAct->getFacetNormal();
1377	for(int ii=0; ii<this->numVars; ii++)
1378	{
1379	dd_set_si(ddineq->matrix[0][ii+1],(*comp)[ii]);
1380	dd_set_si(ddineq->matrix[1][ii+1],1); //Assure we search in the pos. orthant
1381	}
1382	set_addelem(ddineq->linset,1); //We want equality in the first row
1383	//dd_WriteMatrix(stdout,ddineq);
1384	interiorPoint(ddineq,*comp);
1385	/**/
1386	#ifdef gfan_DEBUG
1387	cout << "IP is";
1388	comp->show(); cout << endl;
1389	#endif
1390	//Store the interior point and the UCN
1391	fListPtr->setInteriorPoint( comp );
1392	fListPtr->setUCN( gc.getUCN() );
1393
1394	dd_FreeMatrix(ddineq);
1395	fAct=fAct->next; //iterate
1396	}
1397	}
1398	else/Facets of facets/
1399	{
1400	dd_MatrixPtr ddineq,P,ddakt;
1401	dd_rowset impl_linset, redset;
1402	dd_ErrorType err;
1403	dd_rowindex newpos;
1404
1405	fAct = fListPtr;
1406
1407	#ifdef gfan_DEBUG
1408	cout << "before f2M" << endl;
1409	gc.showFacets();
1410	ddineq = facets2Matrix(gc);
1411	cout << "after f2M" << endl;
1412	gc.showFacets();
1413	#endif
1414
1415	/Now set appropriate linearity/
1416	dd_PolyhedraPtr ddpolyh;
1417	for (int ii=0; ii<this->numFacets; ii++)
1418	{
1419	cout << "------------" << endl;
1420	ddakt = dd_CopyMatrix(ddineq);
1421	set_addelem(ddakt->linset,ii+1);
1422
1423	dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);
1424
1425	dd_WriteMatrix(stdout,ddakt);
1426	ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
1427	P=dd_CopyGenerators(ddpolyh);
1428	dd_WriteMatrix(stdout,P);
1429
1430	/* We loop through each row of P
1431	* normalize it by making all entries integer ones
1432	* and add the resulting vector to the int matrix facet::codim2Facets
1433	*/
1434	for (int jj=1;jj<=P->rowsize;jj++)
1435	{
1436	intvec *n = new intvec(this->numVars);
1437	normalize(P,jj,*n);
1438	//fAct->addCodim2Facet(n);
1439	n->show();
1440	delete n;
1441	}
1442
1443	//intvec *load = new intvec(this->numVars);
1444
1445	dd_FreeMatrix(ddakt);
1446	dd_FreePolyhedra(ddpolyh);
1447	}
1448	//gc.showFacets();
1449
1450	}
1451
1452
1453	//NOTE Hm, comment in and get a crash for free...
1454	//dd_free_global_constants();
1455	gc.writeConeToFile(gc);
1456
1457	/*2nd step
1458	Choose a facet from fListPtr, flip it and forget the previous cone
1459	*/
1460	fAct = fListPtr;
1461	//gcone *gcTmp = new gcone(gc); //copy
1462	//gcTmp->flip(gcTmp->gcBasis,fAct);
1463	//NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
1464
1465	}//void noRevS(gcone &gc)
1466
1467	void normalize(dd_MatrixPtr const &M, int line, intvec &n)
1468	{
1469	mpz_t denom[this->numVars];
1470	for(int ii=0;ii<this->numVars;ii++)
1471	{
1472	mpz_init_set_str(denom[ii],"0",10);
1473	}
1474
1475	mpz_t kgV,tmp;
1476	mpz_init(kgV);
1477	mpz_init(tmp);
1478	//intvec *ivres = new intvec(this->numVars);
1479	// intvec ivres(this->numVars);
1480
1481	for (int ii=0;ii<(M->colsize)-1;ii++)
1482	{
1483	mpz_t z;
1484	mpz_init(z);
1485	mpq_get_den(z,M->matrix[line-1][ii+1]);
1486	//mpz_out_str(stdout,10,z);
1487	mpz_set( denom[ii], z);
1488	mpz_clear(z);
1489	}
1490	/Compute lcm of the denominators/
1491	mpz_set(tmp,denom[0]);
1492	for (int ii=0;ii<(M->colsize)-1;ii++)
1493	{
1494	mpz_lcm(kgV,tmp,denom[ii]);
1495	}
1496	/Multiply the nominators by kgV/
1497	mpq_t qkgV,res;
1498	mpq_init(qkgV);
1499	// mpq_canonicalize(qkgV);
1500	mpq_set_str(qkgV,"1",10);
1501	// mpq_canonicalize(qkgV);
1502	mpq_init(res);
1503	mpq_set_str(res,"1",10);
1504	// mpq_canonicalize(res);
1505	//mpq_set_z(qkgV,kgV);
1506	mpq_set_num(qkgV,kgV);
1507	//mpq_set_den(qkgV,1);
1508	// mpq_canonicalize(qkgV);
1509	for (int ii=0;ii<(M->colsize)-1;ii++)
1510	{
1511	mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
1512	// (*ivres)[ii]=(int)mpz_get_d(mpq_numref(res));
1513	n[ii]=(int)mpz_get_d(mpq_numref(res));
1514	}
1515	//mpz_clear(denom[this->numVars]);
1516	mpz_clear(kgV);
1517	mpq_clear(qkgV); mpq_clear(res);
1518
1519	//return ivres;
1520	}
1521
1522	/** \brief Construct a dd_MatrixPtr from a cone's list of facets
1523	*
1524	*/
1525	dd_MatrixPtr facets2Matrix(gcone const &gc)
1526	{
1527	facet *fAct = new facet();
1528	fAct = gc.facetPtr;
1529
1530	dd_MatrixPtr M;
1531	dd_rowrange ddrows;
1532	dd_colrange ddcols;
1533	ddcols=(this->numVars)+1;
1534	ddrows=this->numFacets;
1535	dd_NumberType numb = dd_Integer;
1536	M=dd_CreateMatrix(ddrows,ddcols);
1537
1538	//dd_set_global_constants();
1539
1540	intvec *comp = new intvec(this->numVars);
1541
1542	for (int jj=0; jj<M->rowsize; jj++)
1543	{
1544	comp = fAct->getFacetNormal();
1545	for (int ii=0; ii<this->numVars; ii++)
1546	{
1547	dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
1548	}
1549	fAct = fAct->next;
1550	}//jj
1551
1552	//delete fAct;
1553	//delete comp;
1554	return M;
1555	}
1556
1557	/** \brief Write information about a cone into a file on disk
1558	*
1559	* This methods writes the information needed for the "second" method into a file.
1560	* The file's is divided in sections containing information on
1561	* 1) the ring
1562	* 2) the cone's Groebner Basis
1563	* 3) the cone's facets
1564	* Each line contains exactly one date
1565	* Each section starts with its name in CAPITALS
1566	*/
1567	void writeConeToFile(gcone const &gc, bool usingIntPoints=FALSE)
1568	{
1569	ofstream gcOutputFile("/tmp/cone1.gc");
1570	facet *fAct = new facet();
1571	fAct = gc.facetPtr;
1572
1573	if (!gcOutputFile)
1574	{
1575	cout << "Error opening file for writing in writeConeToFile" << endl;
1576	}
1577	else
1578	{ /*gcOutputFile << "UCN" << endl;
1579	gcOutputFile << this->UCN << endl;*/
1580	gcOutputFile << "RING" << endl;
1581	//Write this->gcBasis into file
1582	gcOutputFile << "GCBASIS" << endl;
1583	for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
1584	{
1585	gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
1586	}
1587
1588	gcOutputFile << "FACETS" << endl;
1589	while(fAct->next!=NULL)
1590	{
1591	intvec *iv = new intvec(gc.numVars);
1592	iv=fAct->getFacetNormal();
1593	for (int ii=0;ii<iv->length();ii++)
1594	{
1595	if (ii<iv->length()-1)
1596	{
1597	gcOutputFile << (*iv)[ii] << ",";
1598	}
1599	else
1600	{
1601	gcOutputFile << (*iv)[ii] << endl;
1602	}
1603	}
1604	fAct=fAct->next;
1605	//delete iv; iv=NULL;
1606	}
1607
1608	gcOutputFile.close();
1609	//delete fAct; fAct=NULL;
1610	}
1611
1612	}//writeConeToFile(gcone const &gc)
1613
1614	/** \brief Reads a cone from a file identified by its number
1615	*/
1616	void readConeFromFile(int gcNum)
1617	{
1618	}
1619
1620	friend class facet;
1621	};//class gcone
1622
1623	ideal gfan(ideal inputIdeal)
1624	{
1625	int numvar = pVariables;
1626
1627	enum searchMethod {
1628	reverseSearch,
1629	noRevS
1630	};
1631
1632	searchMethod method;
1633	method = noRevS;
1634	// method = reverseSearch;
1635
1636	#ifdef gfan_DEBUG
1637	cout << "Now in subroutine gfan" << endl;
1638	#endif
1639	ring inputRing=currRing; // The ring the user entered
1640	ring rootRing; // The ring associated to the target ordering
1641	ideal res;
1642	facet *fRoot;
1643
1644	if (method==reverseSearch)
1645	{
1646
1647	/* Construct a new ring which will serve as our root*/
1648	rootRing=rCopy0(currRing);
1649	rootRing->order[0]=ringorder_lp;
1650
1651	rComplete(rootRing);
1652	rChangeCurrRing(rootRing);
1653
1654	/* Fetch the inputIdeal into our rootRing */
1655	map theMap=(map)idMaxIdeal(1); //evil hack!
1656	theMap->preimage=NULL; //neccessary?
1657	ideal rootIdeal;
1658	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
1659	#ifdef gfan_DEBUG
1660	cout << "Root ideal is " << endl;
1661	idShow(rootIdeal);
1662	cout << "The root ring is " << endl;
1663	rWrite(rootRing);
1664	cout << endl;
1665	#endif
1666
1667	//gcone *gcRoot = new gcone(); //Instantiate the sink
1668	gcone *gcRoot = new gcone(rootRing,rootIdeal);
1669	gcone *gcAct;
1670	gcAct = gcRoot;
1671	gcAct->numVars=pVariables;
1672	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
1673	idShow(gcAct->gcBasis);
1674	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
1675	//gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1676	/*Now it is time to compute the search facets, respectively start the reverse search.
1677	But since we are in the root all facets should be search facets. IS THIS TRUE?
1678	NOTE: Check for flippability is not very sophisticated
1679	*/
1680	//gcAct->reverseSearch(gcAct);
1681	rChangeCurrRing(rootRing);
1682	res=gcRoot->gcBasis;
1683	}//if method==reverSearch
1684
1685	if(method==noRevS)
1686	{
1687	gcone *gcRoot = new gcone(currRing,inputIdeal);
1688	gcone *gcAct;
1689	gcAct = gcRoot;
1690	gcAct->numVars=pVariables;
1691	gcAct->getGB(inputIdeal);
1692	gcAct->getConeNormals(gcAct->gcBasis);
1693	cout << "ding" << endl;
1694	gcAct->noRevS(*gcAct);
1695	res=gcAct->gcBasis;
1696	}
1697
1698	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
1699	The return type will then be of type LIST_CMD
1700	Assume gfan has finished, thus we have enumerated all the cones
1701	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
1702	Groebner Basis and merge this somehow with LIST_CMD
1703	=> Count the cones!
1704	*/
1705	//rChangeCurrRing(rootRing);
1706	//res=gcAct->gcBasis;
1707	//res=gcRoot->gcBasis;
1708	return res;
1709	//return GBlist;
1710	}
1711	/*
1712	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
1713	*/
1714	#endif

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

Download in other formats: