source: git/Singular/fglm.cc @ 7fee876

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