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{
40 FglmOk,
41 FglmHasOne,
42 FglmNoIdeal,
43 FglmNotReduced,
44 FglmNotZeroDim,
45 FglmIncompatibleRings,
46 // for fglmquot:
47 FglmPolyIsOne,
48 FglmPolyIsZero
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
fglmUpdatesource(const ideal sourceIdeal)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
fglmUpdateresult(ideal & result)89 fglmUpdateresult( ideal & result )
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 sring and dring 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 sring to dring.
118 // if the rings are compatible, it returns FglmOk.
119 // Should be called with currRing==sring
120 FglmState
fglmConsistency(ring sring,ring dring,int * vperm)121 fglmConsistency( ring sring, ring dring, int * vperm )
122 {
123 int k;
124 FglmState state= FglmOk;
125
126 if ( rChar(sring) != rChar(dring) )
127 {
128 WerrorS( "rings must have same characteristic" );
129 state= FglmIncompatibleRings;
130 }
131 if ( (sring->OrdSgn != 1) || (dring->OrdSgn != 1) )
132 {
133 WerrorS( "only works for global orderings" );
134 state= FglmIncompatibleRings;
135 }
136 if ( sring->N != dring->N )
137 {
138 WerrorS( "rings must have same number of variables" );
139 state= FglmIncompatibleRings;
140 }
141 if ( rPar(sring) != rPar(dring) )
142 {
143 WerrorS( "rings must have same number of parameters" );
144 state= FglmIncompatibleRings;
145 }
146 if ( state != FglmOk ) return state;
147 // now the rings have the same number of variables resp. parameters.
148 // check if the names of the variables resp. parameters do agree:
149 int nvar = sring->N;
150 int npar = rPar(sring);
151 int * pperm;
152 if ( npar > 0 )
153 pperm= (int *)omAlloc0( (npar+1)*sizeof( int ) );
154 else
155 pperm= NULL;
156 maFindPerm( sring->names, nvar, rParameter(sring), npar,
157 dring->names, nvar, rParameter(dring), npar, vperm, pperm,
158 dring->cf->type);
159 for ( k= nvar; (k > 0) && (state == FglmOk); k-- )
160 if ( vperm[k] <= 0 )
161 {
162 WerrorS( "variable names do not agree" );
163 state= FglmIncompatibleRings;
164 }
165 for ( k= npar-1; (k >= 0) && (state == FglmOk); k-- )
166 if ( pperm[k] >= 0 )
167 {
168 WerrorS( "parameter names do not agree" );
169 state= FglmIncompatibleRings;
170 }
171 if (pperm != NULL) // OB: ????
172 omFreeSize( (ADDRESS)pperm, (npar+1)*sizeof( int ) );
173 if ( state != FglmOk ) return state;
174 // check if both rings are qrings or not
175 if ( sring->qideal != NULL )
176 {
177 if ( dring->qideal == NULL )
178 {
179 WerrorS( "source ring is a qring, destination ring not" );
180 return FglmIncompatibleRings;
181 }
182 // both rings are qrings, now check if both quotients define the same ideal.
183 // check if sring->qideal is contained in dring->qideal:
184 rChangeCurrRing( dring );
185 nMapFunc nMap=n_SetMap(dring->cf, sring->cf );
186 ideal sqind = idInit( IDELEMS( sring->qideal ), 1 );
187 for ( k= IDELEMS( sring->qideal )-1; k >= 0; k-- )
188 (sqind->m)[k]= p_PermPoly( (sring->qideal->m)[k], vperm, sring,
189 dring, nMap);
190 ideal sqindred = kNF( dring->qideal, NULL, sqind );
191 if ( ! idIs0( sqindred ) )
192 {
193 WerrorS( "the quotients do not agree" );
194 state= FglmIncompatibleRings;
195 }
196 idDelete( & sqind );
197 idDelete( & sqindred );
198 rChangeCurrRing( sring );
199 if ( state != FglmOk ) return state;
200 // check if dring->qideal is contained in sring->qideal:
201 int * dsvperm = (int *)omAlloc0( (nvar+1)*sizeof( int ) );
202 maFindPerm( dring->names, nvar, NULL, 0, sring->names, nvar, NULL, 0,
203 dsvperm, NULL, sring->cf->type);
204 nMap=n_SetMap(currRing->cf, dring->cf);
205 ideal dqins = idInit( IDELEMS( dring->qideal ), 1 );
206 for ( k= IDELEMS( dring->qideal )-1; k >= 0; k-- )
207 (dqins->m)[k]=p_PermPoly( (dring->qideal->m)[k], dsvperm, sring,
208 currRing, nMap);
209 ideal dqinsred = kNF( sring->qideal, NULL, dqins );
210 if ( ! idIs0( dqinsred ) )
211 {
212 WerrorS( "the quotients do not agree" );
213 state= FglmIncompatibleRings;
214 }
215 idDelete( & dqins );
216 idDelete( & dqinsred );
217 omFreeSize( (ADDRESS)dsvperm, (nvar+1)*sizeof( int ) );
218 if ( state != FglmOk ) return state;
219 }
220 else
221 {
222 if ( dring->qideal != NULL )
223 {
224 WerrorS( "source ring is a qring, destination ring not" );
225 return FglmIncompatibleRings;
226 }
227 }
228 return FglmOk;
229 }
230
231 // Checks if the ideal "theIdeal" is zero-dimensional and minimal. It does
232 // not check, if it is reduced.
233 // returns FglmOk if we can use theIdeal for CalculateFunctionals (this
234 // function reports an error if theIdeal is not reduced,
235 // so this need not to be tested here)
236 // FglmNotReduced if theIdeal is not minimal
237 // FglmNotZeroDim if it is not zero-dimensional
238 // FglmHasOne if 1 belongs to theIdeal
239 FglmState
fglmIdealcheck(const ideal theIdeal)240 fglmIdealcheck( const ideal theIdeal )
241 {
242 FglmState state = FglmOk;
243 int power;
244 int k;
245 BOOLEAN * purePowers = (BOOLEAN *)omAlloc0( currRing->N*sizeof( BOOLEAN ) );
246
247 for ( k= IDELEMS( theIdeal ) - 1; (state == FglmOk) && (k >= 0); k-- )
248 {
249 poly p = (theIdeal->m)[k];
250 if (p!=NULL)
251 {
252 if( pIsConstant( p ) ) state= FglmHasOne;
253 else if ( (power= pIsPurePower( p )) > 0 )
254 {
255 fglmASSERT( 0 < power && power <= currRing->N, "illegal power" );
256 if ( purePowers[power-1] == TRUE ) state= FglmNotReduced;
257 else purePowers[power-1]= TRUE;
258 }
259 for ( int l = IDELEMS( theIdeal ) - 1; state == FglmOk && l >= 0; l-- )
260 if ( (k != l) && pDivisibleBy( p, (theIdeal->m)[l] ) )
261 state= FglmNotReduced;
262 }
263 }
264 if ( state == FglmOk )
265 {
266 for ( k= currRing->N-1 ; (state == FglmOk) && (k >= 0); k-- )
267 if ( purePowers[k] == FALSE ) state= FglmNotZeroDim;
268 }
269 omFreeSize( (ADDRESS)purePowers, currRing->N*sizeof( BOOLEAN ) );
270 return state;
271 }
272
273 // The main function for the fglm-Algorithm.
274 // Checks the input-data, and calls fglmzero (see fglmzero.cc).
275 // Returns the new groebnerbasis or 0 if an error occoured.
276 BOOLEAN
fglmProc(leftv result,leftv first,leftv second)277 fglmProc( leftv result, leftv first, leftv second )
278 {
279 FglmState state = FglmOk;
280
281 ring destRing = currRing;
282 // ring destRing = currRing;
283 ideal destIdeal = NULL;
284 ring sourceRing = (ring)first->Data();
285 rChangeCurrRing( sourceRing );
286 // ring sourceRing = currRing;
287
288 int * vperm = (int *)omAlloc0( (sourceRing->N+1)*sizeof( int ) );
289 state= fglmConsistency( sourceRing, destRing, vperm );
290 omFreeSize( (ADDRESS)vperm, (sourceRing->N+1)*sizeof(int) );
291
292 if ( state == FglmOk )
293 {
294 idhdl ih = sourceRing->idroot->get( second->Name(), myynest );
295 if ( (ih != NULL) && (IDTYP(ih)==IDEAL_CMD) )
296 {
297 ideal sourceIdeal;
298 if ( sourceRing->qideal != NULL )
299 sourceIdeal= fglmUpdatesource( IDIDEAL( ih ) );
300 else
301 sourceIdeal = IDIDEAL( ih );
302 state= fglmIdealcheck( sourceIdeal );
303 if ( state == FglmOk )
304 {
305 // Now the settings are compatible with FGLM
306 assumeStdFlag( (leftv)ih );
307 if ( fglmzero( sourceRing, sourceIdeal, destRing, destIdeal, FALSE, (currRing->qideal != NULL) ) == FALSE )
308 state= FglmNotReduced;
309 }
310 } else state= FglmNoIdeal;
311 }
312 if ( currRing != destRing )
313 rChangeCurrRing( destRing );
314 switch (state)
315 {
316 case FglmOk:
317 if ( currRing->qideal != NULL ) fglmUpdateresult( destIdeal );
318 break;
319 case FglmHasOne:
320 destIdeal= idInit(1,1);
321 (destIdeal->m)[0]= pOne();
322 state= FglmOk;
323 break;
324 case FglmIncompatibleRings:
325 WerrorS( "source ring and current ring are incompatible" );
326 destIdeal= NULL;
327 break;
328 case FglmNoIdeal:
329 Werror( "Can't find ideal %s in source ring", second->Name() );
330 destIdeal= NULL;
331 break;
332 case FglmNotZeroDim:
333 Werror( "The ideal %s has to be 0-dimensional", second->Name() );
334 destIdeal= NULL;
335 break;
336 case FglmNotReduced:
337 Werror( "The ideal %s has to be given by a reduced SB", second->Name() );
338 destIdeal= NULL;
339 break;
340 default:
341 destIdeal= idInit(1,1);
342 }
343
344 result->rtyp = IDEAL_CMD;
345 result->data= (void *)destIdeal;
346 setFlag( result, FLAG_STD );
347 return (state != FglmOk);
348 }
349
fglmQuot(ideal first,poly second)350 ideal fglmQuot( ideal first, poly second )
351 {
352 FglmState state = FglmOk;
353
354 ideal sourceIdeal = first;
355 poly quot = second;
356 ideal destIdeal = NULL;
357
358 state = fglmIdealcheck( sourceIdeal );
359 if ( state == FglmOk )
360 {
361 if ( quot == NULL ) state= FglmPolyIsZero;
362 else if ( pIsConstant( quot ) ) state= FglmPolyIsOne;
363 }
364
365 if ( state == FglmOk )
366 {
367 if ( fglmquot( sourceIdeal, quot, destIdeal ) == FALSE )
368 state= FglmNotReduced;
369 }
370
371 switch (state)
372 {
373 case FglmOk:
374 break;
375 case FglmHasOne:
376 destIdeal= idInit(1,1);
377 (destIdeal->m)[0]= pOne();
378 state= FglmOk;
379 break;
380 case FglmNotZeroDim:
381 WerrorS( "The ideal has to be 0-dimensional" );
382 destIdeal= idInit(1,1);
383 break;
384 case FglmNotReduced:
385 WerrorS( "The poly has to be reduced" );
386 destIdeal= idInit(1,1);
387 break;
388 case FglmPolyIsOne:
389 int k;
390 destIdeal= idInit( IDELEMS(sourceIdeal), 1 );
391 for ( k= IDELEMS( sourceIdeal )-1; k >=0; k-- )
392 (destIdeal->m)[k]= pCopy( (sourceIdeal->m)[k] );
393 state= FglmOk;
394 break;
395 case FglmPolyIsZero:
396 destIdeal= idInit(1,1);
397 (destIdeal->m)[0]= pOne();
398 state= FglmOk;
399 break;
400 default:
401 destIdeal= idInit(1,1);
402 }
403
404 return destIdeal;
405 }
406
407 // fglmQuotProc: Calculate I:f with FGLM methods.
408 // Checks the input-data, and calls fglmquot (see fglmzero.cc).
409 // Returns the new groebnerbasis if I:f or 0 if an error occoured.
410 BOOLEAN
fglmQuotProc(leftv result,leftv first,leftv second)411 fglmQuotProc( leftv result, leftv first, leftv second )
412 {
413 FglmState state = FglmOk;
414
415 // STICKYPROT("quotstart\n");
416 ideal sourceIdeal = (ideal)first->Data();
417 poly quot = (poly)second->Data();
418 ideal destIdeal = NULL;
419
420 state = fglmIdealcheck( sourceIdeal );
421 if ( state == FglmOk )
422 {
423 if ( quot == NULL ) state= FglmPolyIsZero;
424 else if ( pIsConstant( quot ) ) state= FglmPolyIsOne;
425 }
426
427 if ( state == FglmOk )
428 {
429 assumeStdFlag( first );
430 if ( fglmquot( sourceIdeal, quot, destIdeal ) == FALSE )
431 state= FglmNotReduced;
432 }
433
434 switch (state)
435 {
436 case FglmOk:
437 break;
438 case FglmHasOne:
439 destIdeal= idInit(1,1);
440 (destIdeal->m)[0]= pOne();
441 state= FglmOk;
442 break;
443 case FglmNotZeroDim:
444 Werror( "The ideal %s has to be 0-dimensional", first->Name() );
445 destIdeal= NULL;
446 break;
447 case FglmNotReduced:
448 Werror( "The poly %s has to be reduced", second->Name() );
449 destIdeal= NULL;
450 break;
451 case FglmPolyIsOne:
452 int k;
453 destIdeal= idInit( IDELEMS(sourceIdeal), 1 );
454 for ( k= IDELEMS( sourceIdeal )-1; k >=0; k-- )
455 (destIdeal->m)[k]= pCopy( (sourceIdeal->m)[k] );
456 state= FglmOk;
457 break;
458 case FglmPolyIsZero:
459 destIdeal= idInit(1,1);
460 (destIdeal->m)[0]= pOne();
461 state= FglmOk;
462 break;
463 default:
464 destIdeal= idInit(1,1);
465 }
466
467 result->rtyp = IDEAL_CMD;
468 result->data= (void *)destIdeal;
469 setFlag( result, FLAG_STD );
470 // STICKYPROT("quotend\n");
471 return (state != FglmOk);
472 } // fglmQuotProt
473
findUni(ideal first)474 ideal findUni( ideal first )
475 {
476 ideal sourceIdeal;
477 ideal destIdeal = NULL;
478 FglmState state;
479
480 sourceIdeal = first;
481
482 state= fglmIdealcheck( sourceIdeal );
483 if ( state == FglmOk )
484 {
485 // check for special cases: if the input contains
486 // univariate polys, try to reduce the problem
487 int i,k;
488 int count=0;
489 BOOLEAN * purePowers = (BOOLEAN *)omAlloc0( currRing->N*sizeof( BOOLEAN ) );
490 for ( k= IDELEMS( sourceIdeal ) - 1; k >= 0; k-- )
491 {
492 if((i=pIsUnivariate(sourceIdeal->m[k]))>0)
493 {
494 if (purePowers[i-1]==0)
495 {
496 purePowers[i-1]=k;
497 count++;
498 if (count==currRing->N) break;
499 }
500 }
501 }
502 if (count==currRing->N)
503 {
504 destIdeal=idInit(currRing->N,1);
505 for(k=currRing->N-1; k>=0; k--) destIdeal->m[k]=pCopy(sourceIdeal->m[purePowers[k]]);
506 }
507 omFreeSize((ADDRESS)purePowers, currRing->N*sizeof( BOOLEAN ) );
508 if (destIdeal!=NULL)
509 state = FglmOk;
510 else if ( FindUnivariateWrapper( sourceIdeal, destIdeal ) == FALSE )
511 state = FglmNotReduced;
512 }
513 switch (state)
514 {
515 case FglmOk:
516 break;
517 case FglmHasOne:
518 destIdeal= idInit(1,1);
519 (destIdeal->m)[0]= pOne();
520 state= FglmOk;
521 break;
522 case FglmNotZeroDim:
523 WerrorS( "The ideal has to be 0-dimensional" );
524 destIdeal= idInit(1,1);
525 break;
526 case FglmNotReduced:
527 Werror( "The ideal has to be reduced" );
528 destIdeal= idInit(1,1);
529 break;
530 default:
531 destIdeal= idInit(1,1);
532 }
533
534 return destIdeal;
535 }
536 // The main function for finduni().
537 // Checks the input-data, and calls FindUnivariateWrapper (see fglmzero.cc).
538 // Returns an ideal containing the univariate Polynomials or 0 if an error
539 // has occoured.
540 BOOLEAN
findUniProc(leftv result,leftv first)541 findUniProc( leftv result, leftv first )
542 {
543 ideal sourceIdeal;
544 ideal destIdeal = NULL;
545 FglmState state;
546
547 sourceIdeal = (ideal)first->Data();
548
549 assumeStdFlag( first );
550 state= fglmIdealcheck( sourceIdeal );
551 if ( state == FglmOk )
552 {
553 // check for special cases: if the input contains
554 // univariate polys, try to reduce the problem
555 int i,k;
556 int count=0;
557 BOOLEAN * purePowers = (BOOLEAN *)omAlloc0( currRing->N*sizeof( BOOLEAN ) );
558 for ( k= IDELEMS( sourceIdeal ) - 1; k >= 0; k-- )
559 {
560 if((i=pIsUnivariate(sourceIdeal->m[k]))>0)
561 {
562 if (purePowers[i-1]==0)
563 {
564 purePowers[i-1]=k;
565 count++;
566 if (count==currRing->N) break;
567 }
568 }
569 }
570 if (count==currRing->N)
571 {
572 destIdeal=idInit(currRing->N,1);
573 for(k=currRing->N-1; k>=0; k--) destIdeal->m[k]=pCopy(sourceIdeal->m[purePowers[k]]);
574 }
575 omFreeSize((ADDRESS)purePowers, currRing->N*sizeof( BOOLEAN ) );
576 if (destIdeal!=NULL)
577 state = FglmOk;
578 else if ( FindUnivariateWrapper( sourceIdeal, destIdeal ) == FALSE )
579 state = FglmNotReduced;
580 }
581 switch (state)
582 {
583 case FglmOk:
584 break;
585 case FglmHasOne:
586 destIdeal= idInit(1,1);
587 (destIdeal->m)[0]= pOne();
588 state= FglmOk;
589 break;
590 case FglmNotZeroDim:
591 Werror( "The ideal %s has to be 0-dimensional", first->Name() );
592 destIdeal= NULL;
593 break;
594 case FglmNotReduced:
595 Werror( "The ideal %s has to be reduced", first->Name() );
596 destIdeal= NULL;
597 break;
598 default:
599 destIdeal= idInit(1,1);
600 }
601
602 result->rtyp = IDEAL_CMD;
603 result->data= (void *)destIdeal;
604
605 return FALSE;
606 }
607 // ----------------------------------------------------------------------------
608 // Local Variables: ***
609 // compile-command: "make Singular" ***
610 // page-delimiter: "^\\(\\|//!\\)" ***
611 // fold-internal-margins: nil ***
612 // End: ***
613