source: git/Singular/fglm.cc @ 066288

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