source: git/Singular/fglm.cc @ 4508ce5

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