My Project  debian-1:4.1.2-p1+ds-2
fglmhom.cc
Go to the documentation of this file.
1 // emacs edit mode for this file is -*- C++ -*-
2 
3 /****************************************
4 * Computer Algebra System SINGULAR *
5 ****************************************/
6 /*
7 * ABSTRACT - The FGLM-Algorithm extended for homogeneous ideals.
8 * Calculates via the hilbert-function a groebner basis.
9 */
10 
11 
12 
13 
14 
15 #include "kernel/mod2.h"
16 #if 0
17 #include "factoryconf.h"
18 #ifndef NOSTREAMIO
19 #include "iostream.h"
20 #endif
21 // #include "Singular/tok.h"
22 #include "kernel/structs.h"
23 #include "kernel/subexpr.h"
24 #include "kernel/polys.h"
25 #include "kernel/ideals.h"
26 #include "polys/monomials/ring.h"
27 #include "kernel/ipid.h"
28 #include "kernel/ipshell.h"
29 #include "polys/monomials/maps.h"
30 #include "fglm.h"
31 #include "fglmvec.h"
32 #include "fglmgauss.h"
33 #include "misc/intvec.h"
34 #include "kernel/GBEngine/kstd1.h"
37 
38 // obachman: Got rid off those "redefiende messages by includeing fglm.h
39 #include "fglm.h"
40 #if 0
41 #define PROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
42 #define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
43 #define PROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
44 #define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
45 #define fglmASSERT(ignore1,ignore2)
46 #endif
47 
48 struct doublepoly
49 {
50  poly sm;
51  poly dm;
52 };
53 
54 class homogElem
55 {
56 public:
57  doublepoly mon;
58  fglmVector v;
59  fglmVector dv;
60  int basis;
61  int destbasis;
62  BOOLEAN inDest;
63  homogElem() : v(), dv(), basis(0), destbasis(0), inDest(FALSE) {}
64  homogElem( poly m, int b, BOOLEAN ind ) :
65  basis(b), inDest(ind)
66  {
67  mon.dm= m;
68  mon.sm= NULL;
69  }
70 #ifndef HAVE_EXPLICIT_CONSTR
71  void initialize( poly m, int b, BOOLEAN ind )
72  {
73  basis = b;
74  inDest = ind;
75  mon.dm = m;
76  mon.sm = NULL;
77  }
78  void initialize( const homogElem h )
79  {
80  basis = h.basis;
81  inDest = h.inDest;
82  mon.dm = h.mon.dm;
83  mon.sm = h.mon.sm;
84  }
85 #endif
86 };
87 
88 struct homogData
89 {
90  ideal sourceIdeal;
91  doublepoly * sourceHeads;
92  int numSourceHeads;
93  ideal destIdeal;
94  int numDestPolys;
95  homogElem * monlist;
96  int monlistmax;
97  int monlistblock;
98  int numMonoms;
99  int basisSize;
100  int overall; // nur zum testen.
101  int numberofdestbasismonoms;
102 // homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL),
103 // numMonoms(0), basisSize(0) {}
104 };
105 
106 int
107 hfglmNextdegree( intvec * source, ideal current, int & deg )
108 {
109  int numelems;
110  intvec * newhilb = hHstdSeries( current, NULL, currRing->qideal );
111 
112  loop
113  {
114  if ( deg < newhilb->length() )
115  {
116  if ( deg < source->length() )
117  numelems= (*newhilb)[deg]-(*source)[deg];
118  else
119  numelems= (*newhilb)[deg];
120  }
121  else
122  {
123  if (deg < source->length())
124  numelems= -(*source)[deg];
125  else
126  {
127  deg= 0;
128  return 0;
129  }
130  }
131  if (numelems != 0)
132  return numelems;
133  deg++;
134  }
135  delete newhilb;
136 }
137 
138 void
139 generateMonoms( poly m, int var, int deg, homogData * dat )
140 {
141  if ( var == (currRing->N) ) {
142  BOOLEAN inSource = FALSE;
143  BOOLEAN inDest = FALSE;
144  poly mon = pCopy( m );
145  pSetExp( mon, var, deg );
146  pSetm( mon );
147  ++dat->overall;
148  int i;
149  for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) {
150  if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) {
151  inSource= TRUE;
152  }
153  }
154  for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) {
155  if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) {
156  inDest= TRUE;
157  }
158  }
159  if ( (!inSource) || (!inDest) ) {
160  int basis = 0;
161  if ( !inSource )
162  basis= ++(dat->basisSize);
163  if ( !inDest )
164  ++dat->numberofdestbasismonoms;
165  if ( dat->numMonoms == dat->monlistmax ) {
166  int k;
167 #ifdef HAVE_EXPLICIT_CONSTR
168  // Expand array using Singulars ReAlloc function
169  dat->monlist=
170  (homogElem * )omReallocSize( dat->monlist,
171  (dat->monlistmax)*sizeof( homogElem ),
172  (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) );
173  for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
174  dat->monlist[k].homogElem();
175 #else
176  // Expand array by generating new one and copying
177  int newsize = dat->monlistmax + dat->monlistblock;
178  homogElem * tempelem = new homogElem[ newsize ];
179  // Copy old elements
180  for ( k= dat->monlistmax - 1; k >= 0; k-- )
181  tempelem[k].initialize( dat->monlist[k] );
182  delete [] homogElem;
183  homogElem = tempelem;
184 #endif
185  dat->monlistmax+= dat->monlistblock;
186  }
187 #ifdef HAVE_EXPLICIT_CONSTR
188  dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
189 #else
190  dat->monlist[dat->numMonoms].initialize( mon, basis, inDest );
191 #endif
192  dat->numMonoms++;
193  if ( inSource && ! inDest ) PROT( "\\" );
194  if ( ! inSource && inDest ) PROT( "/" );
195  if ( ! inSource && ! inDest ) PROT( "." );
196  }
197  else {
198  pDelete( & mon );
199  }
200  return;
201  }
202  else {
203  poly newm = pCopy( m );
204  while ( deg >= 0 ) {
205  generateMonoms( newm, var+1, deg, dat );
206  pIncrExp( newm, var );
207  pSetm( newm );
208  deg--;
209  }
210  pDelete( & newm );
211  }
212  return;
213 }
214 
215 void
216 mapMonoms( ring oldRing, homogData & dat )
217 {
218  int * vperm = (int *)omAlloc( (currRing->N + 1)*sizeof(int) );
219  maFindPerm( oldRing->names, oldRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL );
220  //nSetMap( oldRing->ch, oldRing->parameter, oldRing->P, oldRing->minpoly );
221  nSetMap( oldRing );
222  int s;
223  for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
224 // dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 );
225  // obachman: changed the folowing to reflect the new calling interface of
226  // pPermPoly -- Tim please check whether this is correct!
227  dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 );
228  }
229 }
230 
231 void
232 getVectorRep( homogData & dat )
233 {
234  // Calculate the NormalForms
235  int s;
236  for ( s= 0; s < dat.numMonoms; s++ ) {
237  if ( dat.monlist[s].inDest == FALSE ) {
238  fglmVector v;
239  if ( dat.monlist[s].basis == 0 ) {
240  v= fglmVector( dat.basisSize );
241  // now the monom is in L(source)
242  PROT( "(" );
243  poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm );
244  PROT( ")" );
245  poly temp = nf;
246  while (temp != NULL ) {
247  int t;
248  for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
249  if ( dat.monlist[t].basis > 0 ) {
250  if ( pLmEqual( dat.monlist[t].mon.sm, temp ) ) {
251  number coeff= nCopy( pGetCoeff( temp ) );
252  v.setelem( dat.monlist[t].basis, coeff );
253  }
254  }
255  }
256  pIter(temp);
257  }
258  pDelete( & nf );
259  }
260  else {
261  PROT( "." );
262  v= fglmVector( dat.basisSize, dat.monlist[s].basis );
263  }
264  dat.monlist[s].v= v;
265  }
266  }
267 }
268 
269 void
270 remapVectors( ring oldring, homogData & dat )
271 {
272  //nSetMap( oldring->ch, oldring->parameter, oldring->P, oldring->minpoly );
273  nSetMap( oldring );
274  int s;
275  for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
276  if ( dat.monlist[s].inDest == FALSE ) {
277  int k;
278  fglmVector newv( dat.basisSize );
279  for ( k= dat.basisSize; k > 0; k-- ){
280  number newnum= nMap( dat.monlist[s].v.getelem( k ) );
281  newv.setelem( k, newnum );
282  }
283  dat.monlist[s].dv= newv;
284  }
285  }
286 }
287 
288 void
289 gaussreduce( homogData & dat, int maxnum, int BS )
290 {
291  int s;
292  int found= 0;
293 
294  int destbasisSize = 0;
295  gaussReducer gauss( dat.basisSize );
296 
297  for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) {
298  if ( dat.monlist[s].inDest == FALSE ) {
299  if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) {
300  destbasisSize++;
301  dat.monlist[s].destbasis= destbasisSize;
302  gauss.store();
303  PROT( "." );
304  }
305  else {
306  fglmVector p= gauss.getDependence();
307  poly result = pCopy( dat.monlist[s].mon.dm );
308  pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) );
309  int l = 0;
310  int k;
311  for ( k= 1; k < p.size(); k++ ) {
312  if ( ! p.elemIsZero( k ) ) {
313  while ( dat.monlist[l].destbasis != k )
314  l++;
315  poly temp = pCopy( dat.monlist[l].mon.dm );
316  pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
317  result= pAdd( result, temp );
318  }
319  }
320  if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
321 // PROT2( "(%s)", pString( result ) );
322  PROT( "+" );
323  found++;
324  (dat.destIdeal->m)[dat.numDestPolys]= result;
325  dat.numDestPolys++;
326  if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
327  pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS );
328  IDELEMS( dat.destIdeal )+= BS;
329  }
330 
331  }
332 
333  }
334  }
335  PROT2( "(%i", s );
336  PROT2( "/%i)", dat.numberofdestbasismonoms );
337 }
338 
339 
340 BOOLEAN
341 fglmhomog( ring sourceRing, ideal sourceIdeal, ring destRing, ideal & destIdeal )
342 {
343 #define groebnerBS 16
344  int numGBelems;
345  int deg = 0;
346 
347  homogData dat;
348 
349  // get the hilbert series and the leading monomials of the sourceIdeal:
350  rChangeCurrRing( sourceRing );
351 
352  intvec * hilb = hHstdSeries( sourceIdeal, NULL, currRing->qideal );
353  int s;
354  dat.sourceIdeal= sourceIdeal;
355  dat.sourceHeads= (doublepoly *)omAlloc( IDELEMS( sourceIdeal ) * sizeof( doublepoly ) );
356  for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- )
357  {
358  dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] );
359  }
360  dat.numSourceHeads= IDELEMS( sourceIdeal );
361  rChangeCurrRing( destRing );
362 
363  // Map the sourceHeads to the destRing
364  int * vperm = (int *)omAlloc( (sourceRing->N + 1)*sizeof(int) );
365  maFindPerm( sourceRing->names, sourceRing->N, NULL, 0, currRing->names,
366  currRing->N, NULL, 0, vperm, NULL, currRing->ch);
367  //nSetMap( sourceRing->ch, sourceRing->parameter, sourceRing->P, sourceRing->minpoly );
368  nSetMap( sourceRing );
369  for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- )
370  {
371  dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 );
372  }
373 
374  dat.destIdeal= idInit( groebnerBS, 1 );
375  dat.numDestPolys= 0;
376 
377  while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 )
378  {
379  int num = 0; // the number of monoms of degree deg
380  PROT2( "deg= %i ", deg );
381  PROT2( "num= %i\ngen>", numGBelems );
382  dat.monlistblock= 512;
383  dat.monlistmax= dat.monlistblock;
384 #ifdef HAVE_EXPLICIT_CONSTR
385  dat.monlist= (homogElem *)omAlloc( dat.monlistmax*sizeof( homogElem ) );
386  int j;
387  for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem();
388 #else
389  dat.monlist = new homogElem[ dat.monlistmax ];
390 #endif
391  dat.numMonoms= 0;
392  dat.basisSize= 0;
393  dat.overall= 0;
394  dat.numberofdestbasismonoms= 0;
395 
396  poly start= pOne();
397  generateMonoms( start, 1, deg, &dat );
398  pDelete( & start );
399 
400  PROT2( "(%i/", dat.basisSize );
401  PROT2( "%i)\nvec>", dat.overall );
402  // switch to sourceRing and map monoms
403  rChangeCurrRing( sourceRing );
404  mapMonoms( destRing, dat );
405  getVectorRep( dat );
406 
407  // switch to destination Ring and remap the vectors
408  rChangeCurrRing( destRing );
409  remapVectors( sourceRing, dat );
410 
411  PROT( "<\nred>" );
412  // now do gaussian reduction
413  gaussreduce( dat, numGBelems, groebnerBS );
414 
415 #ifdef HAVE_EXPLICIT_CONSTR
416  omFreeSize( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) );
417 #else
418  delete [] dat.monlist;
419 #endif
420  PROT( "<\n" );
421  }
422  PROT( "\n" );
423  destIdeal= dat.destIdeal;
424  idSkipZeroes( destIdeal );
425  return TRUE;
426 }
427 
428 /* ideal fglmhomProc(leftv first, leftv second) // TODO: Move to Singular/
429 {
430  idhdl dest= currRingHdl;
431  ideal result;
432  // in den durch das erste Argument angegeben Ring schalten:
433  rSetHdl( (idhdl)first->data, TRUE );
434  idhdl ih= currRing->idroot->get( second->name, myynest );
435  ASSERT( ih!=NULL, "Error: Can't find ideal in ring");
436  rSetHdl( dest, TRUE );
437 
438  ideal i= IDIDEAL(ih);
439  fglmhomog( IDRING((idhdl)first->data), i, IDRING(dest), result );
440 
441  return( result );
442 } */
443 
444 #endif
445 
446 // Questions:
447 // Muss ich einen intvec freigeben?
448 // ----------------------------------------------------------------------------
449 // Local Variables: ***
450 // compile-command: "make Singular" ***
451 // page-delimiter: "^\\( \\|//!\\)" ***
452 // fold-internal-margins: nil ***
453 // End: ***
pDivisibleBy
#define pDivisibleBy(a, b)
returns TRUE, if leading monom of a divides leading monom of b i.e., if there exists a expvector c > ...
Definition: polys.h:132
FALSE
#define FALSE
Definition: auxiliary.h:96
pLmEqual
#define pLmEqual(p1, p2)
Definition: polys.h:107
j
int j
Definition: facHensel.cc:105
k
int k
Definition: cfEzgcd.cc:92
rChangeCurrRing
void rChangeCurrRing(ring r)
Definition: polys.cc:15
result
return result
Definition: facAbsBiFact.cc:76
pEnlargeSet
void pEnlargeSet(poly **p, int l, int increment)
Definition: p_polys.cc:3653
ADDRESS
void * ADDRESS
Definition: auxiliary.h:135
h
STATIC_VAR Poly * h
Definition: janet.cc:971
num
CanonicalForm num(const CanonicalForm &f)
Definition: canonicalform.h:330
polys.h
length
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:263
PROT
#define PROT(msg)
Definition: fglm.h:19
pNeg
#define pNeg(p)
Definition: polys.h:186
pDelete
#define pDelete(p_ptr)
Definition: polys.h:175
fglmVector
Definition: fglmvec.h:17
loop
#define loop
Definition: structs.h:79
b
CanonicalForm b
Definition: cfModGcd.cc:4044
pIncrExp
#define pIncrExp(p, i)
Definition: polys.h:42
nf
Definition: gnumpfl.cc:25
found
bool found
Definition: facFactorize.cc:56
nGreaterZero
#define nGreaterZero(n)
Definition: numbers.h:26
stairc.h
TRUE
#define TRUE
Definition: auxiliary.h:100
i
int i
Definition: cfEzgcd.cc:125
omFreeSize
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:258
BOOLEAN
int BOOLEAN
Definition: auxiliary.h:87
structs.h
idSkipZeroes
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
Definition: simpleideals.cc:180
PROT2
#define PROT2(msg, arg)
Definition: fglm.h:21
CxxTest::initialize
void initialize()
mod2.h
pOne
#define pOne()
Definition: polys.h:299
fglmgauss.h
intvec
Definition: intvec.h:18
gaussReducer
Definition: fglmgauss.h:17
pIter
#define pIter(p)
Definition: monomials.h:34
omAlloc
#define omAlloc(size)
Definition: omAllocDecl.h:208
hHstdSeries
intvec * hHstdSeries(ideal S, intvec *modulweight, intvec *wdegree, ideal Q, ring tailRing)
Definition: hilb.cc:1328
maps.h
intvec.h
kNF
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:2823
pAdd
#define pAdd(p, q)
Definition: polys.h:191
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition: polys.h:30
ring.h
kstd1.h
ftmpl_list.h
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:34
m
int m
Definition: cfEzgcd.cc:121
NULL
#define NULL
Definition: omList.c:11
pSetm
#define pSetm(p)
Definition: polys.h:256
maFindPerm
void maFindPerm(char const *const *const preim_names, int preim_n, char const *const *const preim_par, int preim_p, char const *const *const names, int n, char const *const *const par, int nop, int *perm, int *par_perm, n_coeffType ch)
Definition: maps.cc:162
ideals.h
l
int l
Definition: cfEzgcd.cc:93
pSetExp
#define pSetExp(p, i, v)
Definition: polys.h:41
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
p
int p
Definition: cfModGcd.cc:4019
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
fglmvec.h
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
pCopy
#define pCopy(p)
return a copy of the poly
Definition: polys.h:174
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:23
pHead
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:65
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:41
nCopy
#define nCopy(n)
Definition: numbers.h:14
nSetMap
#define nSetMap(R)
Definition: numbers.h:42
if
if(yy_init)
Definition: libparse.cc:1419
omReallocSize
#define omReallocSize(addr, o_size, size)
Definition: omAllocDecl.h:218