source: git/Singular/fglm.cc @ 341696

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