fglm.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 plus extension
8 * Calculate a reduced groebner basis for one ordering, given a
9 * reduced groebner basis for another ordering.
10 * In this file the input is checked. Furthermore we decide, if
11 * the input is 0-dimensional ( then fglmzero.cc is used ) or
12 * if the input is homogeneous ( then fglmhom.cc is used. Yet
13 * not implemented ).
14 * The extension (finduni) finds minimal univariate Polynomials
15 * lying in a 0-dimensional ideal.
16 */
17 
18 
19 
20 
21 #include <kernel/mod2.h>
22 
23 #include <omalloc/omalloc.h>
24 #include <misc/options.h>
25 
26 #include <polys/monomials/ring.h>
27 #include <polys/monomials/maps.h>
28 
29 #include <kernel/polys.h>
30 #include <kernel/ideals.h>
31 
32 #include <kernel/GBEngine/kstd1.h>
33 #include <kernel/fglm/fglm.h>
34 
35 #include <Singular/fglm.h>
36 #include <Singular/ipid.h>
37 #include <Singular/ipshell.h>
38 #include <Singular/tok.h>
39 
40 
41 // internal Version: 1.18.1.6
42 // enumeration to handle the various errors to occour.
43 enum FglmState{
50  // for fglmquot:
53 };
54 
55 // Has to be called, if currRing->qideal != NULL. ( i.e. qring-case )
56 // Then a new ideal is build, consisting of the generators of sourceIdeal
57 // and the generators of currRing->qideal, which are completely reduced by
58 // the sourceIdeal. This means: If sourceIdeal is reduced, then the new
59 // ideal will be reduced as well.
60 // Assumes that currRing == sourceRing
61 ideal fglmUpdatesource( const ideal sourceIdeal )
62 {
63  int k, l, offset;
64  BOOLEAN found;
65  ideal newSource= idInit( IDELEMS( sourceIdeal ) + IDELEMS( currRing->qideal ), 1 );
66  for ( k= IDELEMS( sourceIdeal )-1; k >=0; k-- )
67  (newSource->m)[k]= pCopy( (sourceIdeal->m)[k] );
68  offset= IDELEMS( sourceIdeal );
69  for ( l= IDELEMS( currRing->qideal )-1; l >= 0; l-- )
70  {
71  if ( (currRing->qideal->m)[l] != NULL )
72  {
73  found= FALSE;
74  for ( k= IDELEMS( sourceIdeal )-1; (k >= 0) && (found == FALSE); k-- )
75  if ( pDivisibleBy( (sourceIdeal->m)[k], (currRing->qideal->m)[l] ) )
76  found= TRUE;
77  if ( ! found )
78  {
79  (newSource->m)[offset]= pCopy( (currRing->qideal->m)[l] );
80  offset++;
81  }
82  }
83  }
84  idSkipZeroes( newSource );
85  return newSource;
86 }
87 
88 // Has to be called, if currRing->qideal != NULL, i.e. in qring-case.
89 // Gets rid of the elements of result which are contained in
90 // currRing->qideal and skips Zeroes.
91 // Assumes that currRing == destRing
92 void
94 {
95  int k, l;
96  BOOLEAN found;
97  for ( k= IDELEMS( result )-1; k >=0; k-- )
98  {
99  if ( (result->m)[k] != NULL )
100  {
101  found= FALSE;
102  for ( l= IDELEMS( currRing->qideal )-1; (l >= 0) && ( found == FALSE ); l-- )
103  if ( pDivisibleBy( (currRing->qideal->m)[l], (result->m)[k] ) )
104  found= TRUE;
105  if ( found ) pDelete( & ((result->m)[k]) );
106  }
107  }
108  idSkipZeroes( result );
109 }
110 
111 // Checks if the two rings sringHdl and dringHdl are compatible enough to
112 // be used for the fglm. This means:
113 // 1) Same Characteristic, 2) globalOrderings in both rings,
114 // 3) Same number of variables, 4) same number of parameters
115 // 5) variables in one ring are permutated variables of the other one
116 // 6) parameters in one ring are permutated parameters of the other one
117 // 7) either both rings are rings or both rings are qrings
118 // 8) if they are qrings, the quotientIdeals of both must coincide.
119 // vperm must be a vector of length pVariables+1, initialized by 0.
120 // If both rings are compatible, it stores the permutation of the
121 // variables if mapped from sringHdl to dringHdl.
122 // if the rings are compatible, it returns FglmOk.
123 // Should be called with currRing= IDRING( sringHdl );
124 FglmState
125 fglmConsistency( idhdl sringHdl, idhdl dringHdl, int * vperm )
126 {
127  int k;
128  FglmState state= FglmOk;
129  ring dring = IDRING( dringHdl );
130  ring sring = IDRING( sringHdl );
131 
132  if ( rChar(sring) != rChar(dring) )
133  {
134  WerrorS( "rings must have same characteristic" );
135  state= FglmIncompatibleRings;
136  }
137  if ( (sring->OrdSgn != 1) || (dring->OrdSgn != 1) )
138  {
139  WerrorS( "only works for global orderings" );
140  state= FglmIncompatibleRings;
141  }
142  if ( sring->N != dring->N )
143  {
144  WerrorS( "rings must have same number of variables" );
145  state= FglmIncompatibleRings;
146  }
147  if ( rPar(sring) != rPar(dring) )
148  {
149  WerrorS( "rings must have same number of parameters" );
150  state= FglmIncompatibleRings;
151  }
152  if ( state != FglmOk ) return state;
153  // now the rings have the same number of variables resp. parameters.
154  // check if the names of the variables resp. parameters do agree:
155  int nvar = sring->N;
156  int npar = rPar(sring);
157  int * pperm;
158  if ( npar > 0 )
159  pperm= (int *)omAlloc0( (npar+1)*sizeof( int ) );
160  else
161  pperm= NULL;
162  maFindPerm( sring->names, nvar, rParameter(sring), npar,
163  dring->names, nvar, rParameter(dring), npar, vperm, pperm,
164  dring->cf->type);
165  for ( k= nvar; (k > 0) && (state == FglmOk); k-- )
166  if ( vperm[k] <= 0 )
167  {
168  WerrorS( "variable names do not agree" );
169  state= FglmIncompatibleRings;
170  }
171  for ( k= npar-1; (k >= 0) && (state == FglmOk); k-- )
172  if ( pperm[k] >= 0 )
173  {
174  WerrorS( "parameter names do not agree" );
175  state= FglmIncompatibleRings;
176  }
177  if (pperm != NULL) // OB: ????
178  omFreeSize( (ADDRESS)pperm, (npar+1)*sizeof( int ) );
179  if ( state != FglmOk ) return state;
180  // check if both rings are qrings or not
181  if ( sring->qideal != NULL )
182  {
183  if ( dring->qideal == NULL )
184  {
185  Werror( "%s is a qring, current ring not", sringHdl->id );
186  return FglmIncompatibleRings;
187  }
188  // both rings are qrings, now check if both quotients define the same ideal.
189  // check if sring->qideal is contained in dring->qideal:
190  rSetHdl( dringHdl );
191  nMapFunc nMap=n_SetMap(currRing->cf, sring->cf );
192  ideal sqind = idInit( IDELEMS( sring->qideal ), 1 );
193  for ( k= IDELEMS( sring->qideal )-1; k >= 0; k-- )
194  (sqind->m)[k]= p_PermPoly( (sring->qideal->m)[k], vperm, sring,
195  currRing, nMap);
196  ideal sqindred = kNF( dring->qideal, NULL, sqind );
197  if ( ! idIs0( sqindred ) )
198  {
199  WerrorS( "the quotients do not agree" );
200  state= FglmIncompatibleRings;
201  }
202  idDelete( & sqind );
203  idDelete( & sqindred );
204  rSetHdl( sringHdl );
205  if ( state != FglmOk ) return state;
206  // check if dring->qideal is contained in sring->qideal:
207  int * dsvperm = (int *)omAlloc0( (nvar+1)*sizeof( int ) );
208  maFindPerm( dring->names, nvar, NULL, 0, sring->names, nvar, NULL, 0,
209  dsvperm, NULL, sring->cf->type);
210  nMap=n_SetMap(currRing->cf, dring->cf);
211  ideal dqins = idInit( IDELEMS( dring->qideal ), 1 );
212  for ( k= IDELEMS( dring->qideal )-1; k >= 0; k-- )
213  (dqins->m)[k]=p_PermPoly( (dring->qideal->m)[k], dsvperm, sring,
214  currRing, nMap);
215  ideal dqinsred = kNF( sring->qideal, NULL, dqins );
216  if ( ! idIs0( dqinsred ) )
217  {
218  WerrorS( "the quotients do not agree" );
219  state= FglmIncompatibleRings;
220  }
221  idDelete( & dqins );
222  idDelete( & dqinsred );
223  omFreeSize( (ADDRESS)dsvperm, (nvar+1)*sizeof( int ) );
224  if ( state != FglmOk ) return state;
225  }
226  else
227  {
228  if ( dring->qideal != NULL )
229  {
230  Werror( "current ring is a qring, %s not", sringHdl->id );
231  return FglmIncompatibleRings;
232  }
233  }
234  return FglmOk;
235 }
236 
237 // Checks if the ideal "theIdeal" is zero-dimensional and minimal. It does
238 // not check, if it is reduced.
239 // returns FglmOk if we can use theIdeal for CalculateFunctionals (this
240 // function reports an error if theIdeal is not reduced,
241 // so this need not to be tested here)
242 // FglmNotReduced if theIdeal is not minimal
243 // FglmNotZeroDim if it is not zero-dimensional
244 // FglmHasOne if 1 belongs to theIdeal
245 FglmState
246 fglmIdealcheck( const ideal theIdeal )
247 {
248  FglmState state = FglmOk;
249  int power;
250  int k;
251  BOOLEAN * purePowers = (BOOLEAN *)omAlloc0( currRing->N*sizeof( BOOLEAN ) );
252 
253  for ( k= IDELEMS( theIdeal ) - 1; (state == FglmOk) && (k >= 0); k-- )
254  {
255  poly p = (theIdeal->m)[k];
256  if (p!=NULL)
257  {
258  if( pIsConstant( p ) ) state= FglmHasOne;
259  else if ( (power= pIsPurePower( p )) > 0 )
260  {
261  fglmASSERT( 0 < power && power <= currRing->N, "illegal power" );
262  if ( purePowers[power-1] == TRUE ) state= FglmNotReduced;
263  else purePowers[power-1]= TRUE;
264  }
265  for ( int l = IDELEMS( theIdeal ) - 1; state == FglmOk && l >= 0; l-- )
266  if ( (k != l) && pDivisibleBy( p, (theIdeal->m)[l] ) )
267  state= FglmNotReduced;
268  }
269  }
270  if ( state == FglmOk )
271  {
272  for ( k= currRing->N-1 ; (state == FglmOk) && (k >= 0); k-- )
273  if ( purePowers[k] == FALSE ) state= FglmNotZeroDim;
274  }
275  omFreeSize( (ADDRESS)purePowers, currRing->N*sizeof( BOOLEAN ) );
276  return state;
277 }
278 
279 // The main function for the fglm-Algorithm.
280 // Checks the input-data, and calls fglmzero (see fglmzero.cc).
281 // Returns the new groebnerbasis or 0 if an error occoured.
282 BOOLEAN
283 fglmProc( leftv result, leftv first, leftv second )
284 {
285  FglmState state = FglmOk;
286 
287  idhdl destRingHdl = currRingHdl;
288  // ring destRing = currRing;
289  ideal destIdeal = NULL;
290  idhdl sourceRingHdl = (idhdl)first->data;
291  rSetHdl( sourceRingHdl );
292  // ring sourceRing = currRing;
293 
294  int * vperm = (int *)omAlloc0( (currRing->N+1)*sizeof( int ) );
295  state= fglmConsistency( sourceRingHdl, destRingHdl, vperm );
296  omFreeSize( (ADDRESS)vperm, (currRing->N+1)*sizeof(int) );
297 
298  if ( state == FglmOk )
299  {
300  idhdl ih = currRing->idroot->get( second->Name(), myynest );
301  if ( (ih != NULL) && (IDTYP(ih)==IDEAL_CMD) )
302  {
303  ideal sourceIdeal;
304  if ( currRing->qideal != NULL )
305  sourceIdeal= fglmUpdatesource( IDIDEAL( ih ) );
306  else
307  sourceIdeal = IDIDEAL( ih );
308  state= fglmIdealcheck( sourceIdeal );
309  if ( state == FglmOk )
310  {
311  // Now the settings are compatible with FGLM
312  assumeStdFlag( (leftv)ih );
313  if ( fglmzero( IDRING(sourceRingHdl), sourceIdeal, IDRING(destRingHdl), destIdeal, FALSE, (currRing->qideal != NULL) ) == FALSE )
314  state= FglmNotReduced;
315  }
316  } else state= FglmNoIdeal;
317  }
318  if ( currRingHdl != destRingHdl )
319  rSetHdl( destRingHdl );
320  switch (state)
321  {
322  case FglmOk:
323  if ( currRing->qideal != NULL ) fglmUpdateresult( destIdeal );
324  break;
325  case FglmHasOne:
326  destIdeal= idInit(1,1);
327  (destIdeal->m)[0]= pOne();
328  state= FglmOk;
329  break;
331  Werror( "ring %s and current ring are incompatible", first->Name() );
332  destIdeal= NULL;
333  break;
334  case FglmNoIdeal:
335  Werror( "Can't find ideal %s in ring %s", second->Name(), first->Name() );
336  destIdeal= NULL;
337  break;
338  case FglmNotZeroDim:
339  Werror( "The ideal %s has to be 0-dimensional", second->Name() );
340  destIdeal= NULL;
341  break;
342  case FglmNotReduced:
343  Werror( "The ideal %s has to be given by a reduced SB", second->Name() );
344  destIdeal= NULL;
345  break;
346  default:
347  destIdeal= idInit(1,1);
348  }
349 
350  result->rtyp = IDEAL_CMD;
351  result->data= (void *)destIdeal;
352  setFlag( result, FLAG_STD );
353  return (state != FglmOk);
354 }
355 
356 // fglmQuotProc: Calculate I:f with FGLM methods.
357 // Checks the input-data, and calls fglmquot (see fglmzero.cc).
358 // Returns the new groebnerbasis if I:f or 0 if an error occoured.
359 BOOLEAN
361 {
362  FglmState state = FglmOk;
363 
364  // STICKYPROT("quotstart\n");
365  ideal sourceIdeal = (ideal)first->Data();
366  poly quot = (poly)second->Data();
367  ideal destIdeal = NULL;
368 
369  state = fglmIdealcheck( sourceIdeal );
370  if ( state == FglmOk )
371  {
372  if ( quot == NULL ) state= FglmPolyIsZero;
373  else if ( pIsConstant( quot ) ) state= FglmPolyIsOne;
374  }
375 
376  if ( state == FglmOk )
377  {
378  assumeStdFlag( first );
379  if ( fglmquot( sourceIdeal, quot, destIdeal ) == FALSE )
380  state= FglmNotReduced;
381  }
382 
383  switch (state)
384  {
385  case FglmOk:
386  break;
387  case FglmHasOne:
388  destIdeal= idInit(1,1);
389  (destIdeal->m)[0]= pOne();
390  state= FglmOk;
391  break;
392  case FglmNotZeroDim:
393  Werror( "The ideal %s has to be 0-dimensional", first->Name() );
394  destIdeal= NULL;
395  break;
396  case FglmNotReduced:
397  Werror( "The poly %s has to be reduced", second->Name() );
398  destIdeal= NULL;
399  break;
400  case FglmPolyIsOne:
401  int k;
402  destIdeal= idInit( IDELEMS(sourceIdeal), 1 );
403  for ( k= IDELEMS( sourceIdeal )-1; k >=0; k-- )
404  (destIdeal->m)[k]= pCopy( (sourceIdeal->m)[k] );
405  state= FglmOk;
406  break;
407  case FglmPolyIsZero:
408  destIdeal= idInit(1,1);
409  (destIdeal->m)[0]= pOne();
410  state= FglmOk;
411  break;
412  default:
413  destIdeal= idInit(1,1);
414  }
415 
416  result->rtyp = IDEAL_CMD;
417  result->data= (void *)destIdeal;
418  setFlag( result, FLAG_STD );
419  // STICKYPROT("quotend\n");
420  return (state != FglmOk);
421 } // fglmQuotProt
422 
423 // The main function for finduni().
424 // Checks the input-data, and calls FindUnivariateWrapper (see fglmzero.cc).
425 // Returns an ideal containing the univariate Polynomials or 0 if an error
426 // has occoured.
427 BOOLEAN
429 {
430  ideal sourceIdeal;
431  ideal destIdeal = NULL;
432  FglmState state;
433 
434  sourceIdeal = (ideal)first->Data();
435 
436  assumeStdFlag( first );
437  state= fglmIdealcheck( sourceIdeal );
438  if ( state == FglmOk )
439  {
440  // check for special cases: if the input contains
441  // univariate polys, try to reduce the problem
442  int i,k;
443  int count=0;
444  BOOLEAN * purePowers = (BOOLEAN *)omAlloc0( currRing->N*sizeof( BOOLEAN ) );
445  for ( k= IDELEMS( sourceIdeal ) - 1; k >= 0; k-- )
446  {
447  if((i=pIsUnivariate(sourceIdeal->m[k]))>0)
448  {
449  if (purePowers[i-1]==0)
450  {
451  purePowers[i-1]=k;
452  count++;
453  if (count==currRing->N) break;
454  }
455  }
456  }
457  if (count==currRing->N)
458  {
459  destIdeal=idInit(currRing->N,1);
460  for(k=currRing->N-1; k>=0; k--) destIdeal->m[k]=pCopy(sourceIdeal->m[purePowers[k]]);
461  }
462  omFreeSize((ADDRESS)purePowers, currRing->N*sizeof( BOOLEAN ) );
463  if (destIdeal!=NULL)
464  state = FglmOk;
465  else if ( FindUnivariateWrapper( sourceIdeal, destIdeal ) == FALSE )
466  state = FglmNotReduced;
467  }
468  switch (state)
469  {
470  case FglmOk:
471  break;
472  case FglmHasOne:
473  destIdeal= idInit(1,1);
474  (destIdeal->m)[0]= pOne();
475  state= FglmOk;
476  break;
477  case FglmNotZeroDim:
478  Werror( "The ideal %s has to be 0-dimensional", first->Name() );
479  destIdeal= NULL;
480  break;
481  case FglmNotReduced:
482  Werror( "The ideal %s has to be reduced", first->Name() );
483  destIdeal= NULL;
484  break;
485  default:
486  destIdeal= idInit(1,1);
487  }
488 
489  result->rtyp = IDEAL_CMD;
490  result->data= (void *)destIdeal;
491 
492  return FALSE;
493 }
494 // ----------------------------------------------------------------------------
495 // Local Variables: ***
496 // compile-command: "make Singular" ***
497 // page-delimiter: "^\\( \\|//!\\)" ***
498 // fold-internal-margins: nil ***
499 // End: ***
int status int void size_t count
Definition: si_signals.h:59
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
#define pIsPurePower(p)
Definition: polys.h:231
Class used for (list of) interpreter objects.
Definition: subexpr.h:82
#define fglmASSERT(ignore1, ignore2)
Definition: fglm.h:23
poly kNF(ideal F, ideal Q, poly p, int syzComp, int lazyReduce)
Definition: kstd1.cc:2971
FglmState fglmIdealcheck(const ideal theIdeal)
Definition: fglm.cc:246
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
#define FALSE
Definition: auxiliary.h:94
Compatiblity layer for legacy polynomial operations (over currRing)
return P p
Definition: myNF.cc:203
static int rPar(const ring r)
(r->cf->P)
Definition: ring.h:590
FglmState fglmConsistency(idhdl sringHdl, idhdl dringHdl, int *vperm)
Definition: fglm.cc:125
BOOLEAN fglmquot(ideal sourceIdeal, poly quot, ideal &destIdeal)
Definition: fglmzero.cc:1215
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
int rChar(ring r)
Definition: ring.cc:688
#define TRUE
Definition: auxiliary.h:98
BOOLEAN fglmzero(ring sourceRing, ideal &sourceIdeal, ring destRing, ideal &destideal, BOOLEAN switchBack=TRUE, BOOLEAN deleteIdeal=FALSE)
Definition: fglmzero.cc:1190
#define IDIDEAL(a)
Definition: ipid.h:130
void * ADDRESS
Definition: auxiliary.h:115
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:93
static char const ** rParameter(const ring r)
(r->cf->parameter)
Definition: ring.h:616
const char * Name()
Definition: subexpr.h:120
Definition: idrec.h:34
bool found
Definition: facFactorize.cc:56
void * data
Definition: subexpr.h:88
int myynest
Definition: febase.cc:46
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:10
#define IDTYP(a)
Definition: ipid.h:116
poly p_PermPoly(poly p, const int *perm, const ring oldRing, const ring dst, nMapFunc nMap, const int *par_perm, int OldPar, BOOLEAN use_mult)
Definition: p_polys.cc:3985
ideal fglmUpdatesource(const ideal sourceIdeal)
Definition: fglm.cc:61
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:49
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
idhdl currRingHdl
Definition: ipid.cc:65
#define setFlag(A, F)
Definition: ipid.h:110
void fglmUpdateresult(ideal &result)
Definition: fglm.cc:93
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:221
idrec * idhdl
Definition: ring.h:18
BOOLEAN assumeStdFlag(leftv h)
Definition: subexpr.cc:1474
int i
Definition: cfEzgcd.cc:123
#define pOne()
Definition: polys.h:297
#define IDELEMS(i)
Definition: simpleideals.h:24
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition: coeffs.h:725
void idSkipZeroes(ideal ide)
gives an ideal/module the minimal possible size
BOOLEAN fglmQuotProc(leftv result, leftv first, leftv second)
Definition: fglm.cc:360
#define FLAG_STD
Definition: ipid.h:106
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
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:138
#define IDRING(a)
Definition: ipid.h:124
#define pDelete(p_ptr)
Definition: polys.h:169
void * Data()
Definition: subexpr.cc:1137
const char * id
Definition: idrec.h:39
BOOLEAN fglmProc(leftv result, leftv first, leftv second)
Definition: fglm.cc:283
#define pIsUnivariate(p)
Definition: polys.h:232
BOOLEAN FindUnivariateWrapper(ideal source, ideal &destIdeal)
Definition: fglmzero.cc:1233
polyrec * poly
Definition: hilb.h:10
FglmState
Definition: fglm.cc:43
void rSetHdl(idhdl h)
Definition: ipshell.cc:5038
int offset
Definition: libparse.cc:1091
int BOOLEAN
Definition: auxiliary.h:85
BOOLEAN idIs0(ideal h)
returns true if h is the zero ideal
Definition: fglm.cc:44
void Werror(const char *fmt,...)
Definition: reporter.cc:189
#define omAlloc0(size)
Definition: omAllocDecl.h:211
BOOLEAN findUniProc(leftv result, leftv first)
Definition: fglm.cc:428
return result
Definition: facAbsBiFact.cc:76
int l
Definition: cfEzgcd.cc:94
#define pCopy(p)
return a copy of the poly
Definition: polys.h:168