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