source: git/Singular/fglm.cc @ d828d63

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