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