Context Navigation

source: git/kernel/gfan.cc @ 864fc2

spielwiese

Last change on this file since 864fc2 was 864fc2, checked in by Martin Monerjan, 15 years ago
created gcone::restOfDiv(poly, ideal) to compute f^{G} and gcone::ffG(ideal, ideal) to compute f-f^{G}. However checks and construction of zeropolynomial don't work. This rev will hang in restOfDiv line 488 Substituted fast_map stuff by idrCopyR git-svn-id: file:///usr/local/Singular/svn/trunk@11647 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 20.4 KB

Rev	Line
[b3e45f]	1	/*
[6623f3]	2	Compute the Groebner fan of an ideal
[fc4f9a]	3	$Author: monerjan $
[864fc2]	4	$Date: 2009-04-08 14:02:26 $
	5	$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.31 2009-04-08 14:02:26 monerjan Exp $
	6	$Id: gfan.cc,v 1.31 2009-04-08 14:02:26 monerjan Exp $
[b3e45f]	7	*/
	8
[d3398c]	9	#include "mod2.h"
[34fcf8]	10
	11	#ifdef HAVE_GFAN
	12
[d3398c]	13	#include "kstd1.h"
[864fc2]	14	#include "kutil.h"
[50ab25a]	15	#include "intvec.h"
[2c6535]	16	#include "polys.h"
	17	#include "ideals.h"
[198a339]	18	#include "kmatrix.h"
[90adba8]	19	#include "fast_maps.h" //Mapping of ideals
[798f721]	20	#include "maps.h"
[864fc2]	21	#include "ring.h"
	22	#include "prCopy.h"
[aba53c]	23	#include "iostream.h" //deprecated
[fc4f9a]	24
	25	//Hacks for different working places
[54248e]	26	#define ITWM
[fc4f9a]	27
	28	#ifdef UNI
[c5d8dd]	29	#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
	30	#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
[fc4f9a]	31	#endif
	32
	33	#ifdef HOME
	34	#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
	35	#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
	36	#endif
	37
	38	#ifdef ITWM
[aba53c]	39	#include "/u/slg/monerjan/cddlib/include/setoper.h"
	40	#include "/u/slg/monerjan/cddlib/include/cdd.h"
[fc4f9a]	41	#endif
[d3398c]	42
[b3e45f]	43	#ifndef gfan_DEBUG
	44	#define gfan_DEBUG
	45	#endif
[9f9b142]	46
	47	//#include gcone.h
	48
[fe919c]	49	/**
	50	*\brief Class facet
[53e33d9]	51	* Implements the facet structure as a linked list
[fe919c]	52	*
	53	*/
[8bdaab]	54	class facet
	55	{
	56	private:
[f5b077]	57	/** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
	58	intvec *fNormal;
	59
	60	/** \brief The Groebner basis on the other side of a shared facet
	61	*
	62	* In order not to have to compute the flipped GB twice we store the basis we already get
	63	* when identifying search facets. Thus in the next step of the reverse search we can
	64	* just copy the old cone and update the facet and the gcBasis.
	65	* facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
	66	*/
	67	ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip
	68
	69
[8bdaab]	70	public:
[fe919c]	71	/** The default constructor. Do I need a constructor of type facet(intvec)? */
	72	facet()
[fc4f9a]	73	{
[fe919c]	74	// Pointer to next facet. */
	75	/* Defaults to NULL. This way there is no need to check explicitly */
	76	this->next=NULL;
[fc4f9a]	77	}
	78
[fe919c]	79	/** The default destructor */
	80	~facet(){;}
	81
[54248e]	82	/** Stores the facet normal \param intvec*/
	83	void setFacetNormal(intvec *iv){
[a9fa92]	84	this->fNormal = ivCopy(iv);
[54248e]	85	//return;
	86	}
	87
[a9fa92]	88	/** Hopefully returns the facet normal */
	89	intvec *getFacetNormal()
	90	{
	91	//this->fNormal->show(); cout << endl;
	92	return this->fNormal;
	93	}
	94
[54248e]	95	/** Method to print the facet normal*/
	96	void printNormal()
	97	{
	98	fNormal->show();
[fc4f9a]	99	}
	100
[f5b077]	101	/** Store the flipped GB*/
	102	void setFlipGB(ideal I)
	103	{
	104	this->flipGB=I;
	105	}
	106
	107	/** Print the flipped GB*/
	108	void printFlipGB()
	109	{
	110	idShow(this->flipGB);
	111	}
[798f721]	112
[18765e]	113	bool isFlippable; //flippable facet? Want to have cone->isflippable.facet[i]
	114	bool isIncoming; //Is the facet incoming or outgoing?
[aba53c]	115	facet *next; //Pointer to next facet
[8bdaab]	116	};
	117
[fe919c]	118	/**
[90adba8]	119	*\brief Implements the cone structure
	120	*
	121	* A cone is represented by a linked list of facet normals
	122	* @see facet
[fe919c]	123	*/
[8bdaab]	124	/*class gcone
	125	finally this should become s.th. like gconelib.{h,cc} to provide an API
	126	*/
	127	class gcone
	128	{
[18765e]	129	private:
[2b7bb5d]	130	int numFacets; //#of facets of the cone
[18765e]	131
	132	public:
[90adba8]	133	/** \brief Default constructor.
	134	*
	135	* Initialises this->next=NULL and this->facetPtr=NULL
	136	*/
[54248e]	137	gcone()
	138	{
	139	this->next=NULL;
	140	this->facetPtr=NULL;
	141	}
[fc4f9a]	142	~gcone(); //destructor
[2b7bb5d]	143
[53e33d9]	144	/** Pointer to the first facet */
[2b7bb5d]	145	facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals
[aba53c]	146	poly gcMarkedTerm; //marked terms of the cone's Groebner basis
[2b7bb5d]	147	int numVars; //#of variables in the ring
	148
	149	/** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
	150	ideal gcBasis; //GB of the cone, set by gcone::getGB();
[18765e]	151	gcone next; //Pointer to previous* cone in search tree
[2b7bb5d]	152
[90adba8]	153	/** \brief Compute the normals of the cone
	154	*
	155	* This method computes a representation of the cone in terms of facet normals. It takes an ideal
	156	* as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
	157	* Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
	158	* Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
[2b7bb5d]	159	* each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
	160	* the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
[90adba8]	161	* As a result of this procedure the pointer facetPtr points to the first facet of the cone.
	162	*/
[f5b077]	163	void getConeNormals(const ideal &I)
[03de21]	164	{
[90adba8]	165	#ifdef gfan_DEBUG
	166	cout << "* Computing Inequalities... *" << endl;
	167	#endif
	168	//All variables go here - except ineq matrix and *v, see below
	169	int lengthGB=IDELEMS(I); // # of polys in the groebner basis
	170	int pCompCount; // # of terms in a poly
	171	poly aktpoly;
	172	int numvar = pVariables; // # of variables in a polynomial (or ring?)
	173	int leadexp[numvar]; // dirty hack of exp.vects
	174	int aktexp[numvar];
	175	int cols,rows; // will contain the dimensions of the ineq matrix - deprecated by
	176	dd_rowrange ddrows;
	177	dd_colrange ddcols;
	178	dd_rowset ddredrows; // # of redundant rows in ddineq
	179	dd_rowset ddlinset; // the opposite
	180	dd_rowindex ddnewpos; // all to make dd_Canonicalize happy
	181	dd_NumberType ddnumb=dd_Real; //Number type
	182	dd_ErrorType dderr=dd_NoError; //
	183	// End of var declaration
	184	#ifdef gfan_DEBUG
	185	cout << "The Groebner basis has " << lengthGB << " elements" << endl;
	186	cout << "The current ring has " << numvar << " variables" << endl;
	187	#endif
	188	cols = numvar;
[198a339]	189
[90adba8]	190	//Compute the # inequalities i.e. rows of the matrix
	191	rows=0; //Initialization
	192	for (int ii=0;ii<IDELEMS(I);ii++)
[03de21]	193	{
[90adba8]	194	aktpoly=(poly)I->m[ii];
	195	rows=rows+pLength(aktpoly)-1;
[03de21]	196	}
[90adba8]	197	#ifdef gfan_DEBUG
	198	cout << "rows=" << rows << endl;
	199	cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
	200	#endif
	201	dd_rowrange aktmatrixrow=0; // needed to store the diffs of the expvects in the rows of ddineq
	202	dd_set_global_constants();
	203	ddrows=rows;
	204	ddcols=cols;
	205	dd_MatrixPtr ddineq; //Matrix to store the inequalities
	206	ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
	207
	208	// We loop through each g\in GB and compute the resulting inequalities
	209	for (int i=0; i<IDELEMS(I); i++)
	210	{
	211	aktpoly=(poly)I->m[i]; //get aktpoly as i-th component of I
	212	pCompCount=pLength(aktpoly); //How many terms does aktpoly consist of?
	213	cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
	214
	215	int v=(int )omAlloc((numvar+1)*sizeof(int));
	216	pGetExpV(aktpoly,v); //find the exp.vect in v[1],...,v[n], use pNext(p)
	217
	218	//Store leadexp for aktpoly
	219	for (int kk=0;kk<numvar;kk++)
	220	{
	221	leadexp[kk]=v[kk+1];
	222	//Since we need to know the difference of leadexp with the other expvects we do nothing here
	223	//but compute the diff below
	224	}
	225
	226
	227	while (pNext(aktpoly)!=NULL) //move to next term until NULL
	228	{
	229	aktpoly=pNext(aktpoly);
	230	pSetm(aktpoly); //doesn't seem to help anything
	231	pGetExpV(aktpoly,v);
	232	for (int kk=0;kk<numvar;kk++)
	233	{
	234	aktexp[kk]=v[kk+1];
	235	//ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk]; //dito
	236	dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
	237	}
	238	aktmatrixrow=aktmatrixrow+1;
	239	}//while
	240
	241	} //for
	242
	243	//Maybe add another row to contain the constraints of the standard simplex?
[c5d8dd]	244
[90adba8]	245	#ifdef gfan_DEBUG
	246	cout << "The inequality matrix is" << endl;
	247	dd_WriteMatrix(stdout, ddineq);
	248	#endif
[e2fb7a]	249
[90adba8]	250	// The inequalities are now stored in ddineq
	251	// Next we check for superflous rows
	252	ddredrows = dd_RedundantRows(ddineq, &dderr);
	253	if (dderr!=dd_NoError) // did an error occur?
	254	{
	255	dd_WriteErrorMessages(stderr,dderr); //if so tell us
	256	} else
	257	{
	258	cout << "Redundant rows: ";
	259	set_fwrite(stdout, ddredrows); //otherwise print the redundant rows
	260	}//if dd_Error
	261
	262	//Remove reduntant rows here!
	263	dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
	264	ddrows = ddineq->rowsize; //Size of the matrix with redundancies removed
	265	ddcols = ddineq->colsize;
[aba53c]	266	#ifdef gfan_DEBUG
[90adba8]	267	cout << "Having removed redundancies, the normals now read:" << endl;
	268	dd_WriteMatrix(stdout,ddineq);
	269	cout << "Rows = " << ddrows << endl;
	270	cout << "Cols = " << ddcols << endl;
[aba53c]	271	#endif
[0ad81f]	272
[90adba8]	273	/Write the normals into class facet/
	274	#ifdef gfan_DEBUG
	275	cout << "Creating list of normals" << endl;
	276	#endif
	277	/The pointer fRoot should be the return value of this function*/
[a9fa92]	278	facet *fRoot = new facet(); //instantiate new facet
	279	this->facetPtr = fRoot; //set variable facetPtr of class gcone to first facet
[90adba8]	280	facet *fAct; //instantiate pointer to active facet
[0ad81f]	281	fAct = fRoot; //Seems to do the trick. fRoot and fAct have to point to the same adress!
[90adba8]	282	std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
	283	for (int kk = 0; kk<ddrows; kk++)
[aba53c]	284	{
[90adba8]	285	intvec *load = new intvec(numvar); //intvec to store a single facet normal that will then be stored via setFacetNormal
	286	for (int jj = 1; jj <ddcols; jj++)
	287	{
	288	double *foo;
	289	foo = (double*)ddineq->matrix[kk][jj]; //get entry from actual position
	290	#ifdef gfan_DEBUG
	291	std::cout << "fAct is " << *foo << " at " << fAct << std::endl;
	292	#endif
	293	(load)[jj-1] = (int)foo; //store typecasted entry at pos jj-1 of load
	294	//check for flipability here
	295	if (jj<ddcols) //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
	296	{
[0ad81f]	297	//fAct->next = new facet(); //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor
	298	}
	299	}//for (int jj = 1; jj <ddcols; jj++)
	300	/Quick'n'dirty hack for flippability/
	301	bool isFlippable;
	302	for (int jj = 0; jj<this->numVars; jj++)
	303	{
	304	intvec *ivCanonical = new intvec(this->numVars);
	305	(*ivCanonical)[jj]=1;
	306	if (dotProduct(load,ivCanonical)>=0)
	307	{
	308	isFlippable=FALSE;
[90adba8]	309	}
[0ad81f]	310	else
	311	{
	312	isFlippable=TRUE;
	313	}
	314	delete ivCanonical;
	315	}//for (int jj = 0; jj<this->numVars; jj++)
	316	if (isFlippable==FALSE)
	317	{
	318	cout << "Ignoring facet";
	319	load->show();
	320	//fAct->next=NULL;
	321	}
	322	else
	323	{ /Now load should be full and we can call setFacetNormal/
	324	fAct->setFacetNormal(load);
	325	fAct->next = new facet();
	326	//fAct->printNormal();
	327	fAct=fAct->next; //this should definitely not be called in the above while-loop :D
	328	}//if (isFlippable==FALSE)
[a9fa92]	329	delete load;
[0ad81f]	330	}//for (int kk = 0; kk<ddrows; kk++)
[90adba8]	331	/*
	332	Now we should have a linked list containing the facet normals of those facets that are
	333	-irredundant
	334	-flipable
	335	Adressing is done via *facetPtr
	336	*/
	337
	338	//Clean up but don't delete the return value! (Whatever it will turn out to be)
	339	dd_FreeMatrix(ddineq);
	340	set_free(ddredrows);
	341	free(ddnewpos);
	342	set_free(ddlinset);
	343	dd_free_global_constants();
	344
	345	}//method getConeNormals(ideal I)
	346
[2b7bb5d]	347
[798f721]	348	/** \brief Compute the Groebner Basis on the other side of a shared facet
	349	*
	350	* Implements algorithm 4.3.2 from Anders' thesis.
	351	* As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
	352	* suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$ in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
	353	* is parallel to \f$ leadexp(g) \f$
[5cd07b]	354	* Parallelity is checked using basic linear algebra. See gcone::isParallel.
	355	* Other possibilities includes computing the rank of the matrix consisting of the vectors in question and
[2b7bb5d]	356	* computing an interior point of the facet and taking all terms having the same weight with respect
	357	* to this interior point.
[798f721]	358	*\param ideal, facet
[2b7bb5d]	359	* Input: a marked,reduced Groebner basis and a facet
[798f721]	360	*/
[2b7bb5d]	361	void flip(ideal gb, facet *f) //Compute "the other side"
	362	{
[a9fa92]	363	intvec *fNormal = new intvec(this->numVars); //facet normal, check for parallelity
[5cd07b]	364	fNormal = f->getFacetNormal(); //read this->fNormal;
	365	#ifdef gfan_DEBUG
[0ad81f]	366	cout << "===" << endl;
	367	cout << "running gcone::flip" << endl;
[5cd07b]	368	cout << "fNormal=";
[a9fa92]	369	fNormal->show();
[5cd07b]	370	cout << endl;
	371	#endif
[798f721]	372	/1st step: Compute the initial ideal/
[5cd07b]	373	poly initialFormElement[IDELEMS(gb)]; //array of #polys in GB to store initial form
[f5b077]	374	ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
[2b7bb5d]	375	poly aktpoly;
[5cd07b]	376	intvec *check = new intvec(this->numVars); //array to store the difference of LE and v
[a9fa92]	377
[2b7bb5d]	378	for (int ii=0;ii<IDELEMS(gb);ii++)
	379	{
[a9fa92]	380	aktpoly = (poly)gb->m[ii];
[2b7bb5d]	381	int v=(int )omAlloc((this->numVars+1)*sizeof(int));
[5cd07b]	382	int leadExpV=(int )omAlloc((this->numVars+1)*sizeof(int));
[2b7bb5d]	383	pGetExpV(aktpoly,leadExpV); //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
[5cd07b]	384	initialFormElement[ii]=pHead(aktpoly);
[a9fa92]	385
[2b7bb5d]	386	while(pNext(aktpoly)!=NULL) /loop trough terms and check for parallelity/
	387	{
	388	aktpoly=pNext(aktpoly); //next term
[5cd07b]	389	pSetm(aktpoly);
	390	pGetExpV(aktpoly,v);
	391	/* Convert (int)v into (intvec)check */
	392	for (int jj=0;jj<this->numVars;jj++)
[2b7bb5d]	393	{
[5cd07b]	394	//cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
	395	//cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
	396	(*check)[jj]=v[jj+1]-leadExpV[jj+1];
[a9fa92]	397	}
[f5b077]	398	#ifdef gfan_DEBUG
[5cd07b]	399	cout << "check=";
	400	check->show();
	401	cout << endl;
[f5b077]	402	#endif
[5cd07b]	403	if (isParallel(check,fNormal)) //pass *check when
[a9fa92]	404	{
[5cd07b]	405	cout << "Parallel vector found, adding to initialFormElement" << endl;
	406	initialFormElement[ii] = pAdd(initialFormElement[ii],(poly)pHead(aktpoly));
	407	}
	408	}//while
	409	cout << "Initial Form=";
	410	pWrite(initialFormElement[ii]);
	411	cout << "---" << endl;
[f5b077]	412	/Now initialFormElement must be added to (ideal)initialForm /
	413	initialForm->m[ii]=initialFormElement[ii];
[5cd07b]	414	}//for
[f5b077]	415	f->setFlipGB(initialForm);
	416	#ifdef gfan_DEBUG
	417	cout << "Initial ideal is: " << endl;
[0ad81f]	418	idShow(initialForm);
[f5b077]	419	f->printFlipGB();
	420	cout << "===" << endl;
	421	#endif
	422	delete check;
	423
[0ad81f]	424	/2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight/
	425	/*Substep 2.1
	426	compute $G_{-\alpha}(in_v(I))
	427	see journal p. 66
	428	*/
	429	ring srcRing=currRing;
	430	ring dstRing=rCopy0(srcRing);
	431	dstRing->order[0]=ringorder_a;
	432	//tmpring->order[1]=ringorder_dp;
	433	//tmpring->order[2]=ringorder_C;
	434	dstRing->wvhdl[0] =( int )omAlloc((fNormal->length())sizeof(int)); //found in Singular/ipshell.cc
[f5b077]	435
	436	for (int ii=0;ii<this->numVars;ii++)
[0ad81f]	437	{
	438	dstRing->wvhdl[0][ii]=-(*fNormal)[ii]; //What exactly am I doing here?
	439	//cout << tmpring->wvhdl[0][ii] << endl;
[f5b077]	440	}
[0ad81f]	441	rComplete(dstRing);
	442	rChangeCurrRing(dstRing);
[864fc2]	443	//map theMap=(map)idMaxIdeal(1);
[0ad81f]	444	rWrite(currRing); cout << endl;
	445	ideal ina;
[864fc2]	446	//ina=fast_map(initialForm,srcRing,(ideal)theMap,dstRing);
	447	ina=idrCopyR(initialForm,srcRing);
[0ad81f]	448	cout << "ina=";
	449	idShow(ina); cout << endl;
[f5b077]	450	ideal H;
[0ad81f]	451	H=kStd(ina,NULL,isHomog,NULL); //we know it is homogeneous
	452	idSkipZeroes(H);
	453	cout << "H="; idShow(H); cout << endl;
	454	/*Substep 2.2
	455	do the lifting
	456	*/
	457	rChangeCurrRing(srcRing);
	458	ideal srcRing_H;
	459	ideal srcRing_HH;
	460	//map theMap = (map)idMaxIdeal(1);
[864fc2]	461	//srcRing_H=fast_map(H,dstRing,(ideal)theMap,srcRing);
	462	srcRing_H=idrCopyR(H,dstRing);
	463	idShow(srcRing_H);
	464	srcRing_HH=ffG(srcRing_H,this->gcBasis);
	465	idShow(srcRing_HH);
[0ad81f]	466	/*Substep 2.3
	467	turn the minimal basis into a reduced one
	468	*/
[f5b077]	469
[2b7bb5d]	470	}//void flip(ideal gb, facet *f)
[90adba8]	471
[864fc2]	472	/** \brief Compute the remainder of a polynomial by a given ideal
	473	*
	474	* Compute \f$ f^{\mathcal{G}} \f$
	475	* Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
	476	* However, since we are only interested in the remainder, there is no need to
	477	* compute the factors \f$ a_i \f$
	478	*/
	479	poly restOfDiv(poly const &f, ideal const &I)
	480	{
	481	cout << "Entering restOfDiv" << endl;
	482	poly p=f;
	483	pWrite(p);
	484	//poly r=kCreateZeroPoly(,currRing,currRing); //The 0-polynomial, hopefully
	485	poly r=0;
	486	int ii;
	487	bool divOccured;
	488	cout << "TICK 1" << endl; //Hangs after this line. Zeropoly stuff
	489	while (pHead(p)!=NULL)
	490	{
	491	ii=1;
	492	divOccured=TRUE;
	493	while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
	494	{
	495	cout << "TICK 2" << endl;
	496	if (pLmDivisibleBy(I->m[ii],p))
	497	{
	498	cout << "TICK 3" << endl;
	499	p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
	500	cout << "TICK 4" << endl;
	501	pWrite(p);
	502	pSetm(p);
	503	divOccured=TRUE;
	504	}
	505	else
	506	{
	507	ii += 1;
	508	}//if (pLmDivisibleBy(I->m[ii],p,currRing))
	509	if (divOccured==FALSE)
	510	{
	511	cout << "TICK 5" << endl;
	512	r=pAdd(r,pHead(p));
	513	cout << "TICK 6" << endl;
	514	p=pSub(p,pHead(p));
	515	cout << "TICK 7" << endl;
	516	pWrite(p);
	517	}//if (divOccured==FALSE)
	518	}//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
	519	}//while (p!=0)
	520	return r;
	521	}//poly restOfDiv(poly const &f, ideal const &I)
	522
	523	/** \brief Compute \f$ f-f^{\mathcal{G}} \f$
	524	*/
	525	ideal ffG(ideal const &H, ideal const &G)
	526	{
	527	cout << "Entering ffG" << endl;
	528	int size=IDELEMS(G);
	529	ideal res=idInit(size,1);
	530	for (int ii=0;ii<size;ii++)
	531	{
	532	res->m[ii]=restOfDiv(H->m[ii],G);
	533	}
	534	return res;
	535	}
	536
[90adba8]	537	/** \brief Compute a Groebner Basis
	538	*
[2b7bb5d]	539	* Computes the Groebner basis and stores the result in gcone::gcBasis
[90adba8]	540	*\param ideal
	541	*\return void
	542	*/
[0ad81f]	543	void getGB(ideal const &inputIdeal)
[90adba8]	544	{
	545	ideal gb;
	546	gb=kStd(inputIdeal,NULL,testHomog,NULL);
	547	idSkipZeroes(gb);
	548	this->gcBasis=gb; //write the GB into gcBasis
[0ad81f]	549	}//void getGB
[90adba8]	550
	551	ideal GenGrbWlk(ideal, ideal); //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
[e2fb7a]	552
[90adba8]	553
[5cd07b]	554	/** \brief Compute the dot product of two intvecs
	555	*
	556	*/
[0ad81f]	557	int dotProduct(intvec const &iva, intvec const &ivb)
[5cd07b]	558	{
[0ad81f]	559	//intvec iva=a;
	560	//intvec ivb=b;
[5cd07b]	561	int res=0;
	562	for (int i=0;i<this->numVars;i++)
	563	{
	564	res = res+(iva[i]*ivb[i]);
	565	}
	566	return res;
[0ad81f]	567	}//int dotProduct
[90adba8]	568
[5cd07b]	569	/** \brief Check whether two intvecs are parallel
	570	*
	571	* \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
	572	*/
[0ad81f]	573	bool isParallel(intvec const &a, intvec const &b)
[5cd07b]	574	{
	575	int lhs,rhs;
[0ad81f]	576	lhs=dotProduct(a,b)*dotProduct(a,b);
	577	rhs=dotProduct(a,a)*dotProduct(b,b);
[5cd07b]	578	cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
	579	if (lhs==rhs)
	580	{
	581	return TRUE;
	582	}
	583	else
	584	{
	585	return FALSE;
	586	}
	587	}//bool isParallel
	588	};//class gcone
[2c6535]	589
	590	ideal gfan(ideal inputIdeal)
	591	{
[54248e]	592	int numvar = pVariables;
	593
[b3e45f]	594	#ifdef gfan_DEBUG
[fc4f9a]	595	cout << "Now in subroutine gfan" << endl;
[b3e45f]	596	#endif
[90adba8]	597	ring inputRing=currRing; // The ring the user entered
	598	ring rootRing; // The ring associated to the target ordering
[b3e45f]	599	ideal res;
[90adba8]	600	//matrix ineq; //Matrix containing the boundary inequalities
[9f9b142]	601	facet *fRoot;
	602
[90adba8]	603	/*
	604	1. Select target order, say dp.
	605	2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
	606	3. getConeNormals
	607	*/
[9f9b142]	608
[90adba8]	609	/* Construct a new ring which will serve as our root
	610	Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
[864fc2]	611	resolved 07.04.2009 MM
[90adba8]	612	*/
[9f9b142]	613	rootRing=rCopy0(currRing);
[f5b077]	614	rootRing->order[0]=ringorder_lp;
[9f9b142]	615	rComplete(rootRing);
[f5b077]	616	rChangeCurrRing(rootRing);
	617
[90adba8]	618	/* Fetch the inputIdeal into our rootRing */
[f5b077]	619	map theMap=(map)idMaxIdeal(1); //evil hack!
[864fc2]	620	theMap->preimage=NULL; //neccessary?
[f5b077]	621	ideal rootIdeal;
	622	rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
[90adba8]	623	#ifdef gfan_DEBUG
	624	cout << "Root ideal is " << endl;
[f5b077]	625	idShow(rootIdeal);
	626	cout << "The root ring is " << endl;
[9f9b142]	627	rWrite(rootRing);
[90adba8]	628	cout << endl;
	629	#endif
	630
[9f9b142]	631	gcone *gcRoot = new gcone(); //Instantiate the sink
	632	gcone *gcAct;
	633	gcAct = gcRoot;
[2b7bb5d]	634	gcAct->numVars=pVariables;
[f5b077]	635	gcAct->getGB(rootIdeal); //sets gcone::gcBasis
	636	idShow(gcAct->gcBasis);
[90adba8]	637	gcAct->getConeNormals(gcAct->gcBasis); //hopefully compute the normals
[a9fa92]	638	gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
[90adba8]	639	/*Now it is time to compute the search facets, respectively start the reverse search.
[798f721]	640	But since we are in the root all facets should be search facets. IS THIS TRUE?
[90adba8]	641	MIND: AS OF NOW, THE LIST OF FACETS IS NOT PURGED OF NON-FLIPPAPLE FACETS
	642	*/
[9f9b142]	643
[90adba8]	644	/*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
	645	The return type will then be of type LIST_CMD
	646	Assume gfan has finished, thus we have enumerated all the cones
	647	Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
	648	Groebner Basis and merge this somehow with LIST_CMD
	649	=> Count the cones!
[fc4f9a]	650	*/
[0ad81f]	651	rChangeCurrRing(inputRing);
[90adba8]	652	res=gcAct->gcBasis;
	653	//cout << fRoot << endl;
[b3e45f]	654	return res;
[90adba8]	655	//return GBlist;
[d3398c]	656	}
[9f9b142]	657	/*
	658	Since gfan.cc is #included from extra.cc there must not be a int main(){} here
	659	*/
[34fcf8]	660	#endif

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

Download in other formats: