source: git/Singular/fglm.cc @ dcd92d

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