source: git/Singular/fglm.cc @ bee06d

spielwiese
Last change on this file since bee06d was b1dfaf, checked in by Frank Seelisch <seelisch@…>, 14 years ago
patch from Kai (checked for problems under Windows OS: no problems) git-svn-id: file:///usr/local/Singular/svn/trunk@13210 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.2 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 <kernel/mod2.h>
20
21#ifdef HAVE_FACTORY
22#include <Singular/tok.h>
23#include <kernel/options.h>
24#include <kernel/polys.h>
25#include <kernel/ideals.h>
26#include <kernel/ring.h>
27#include <Singular/ipid.h>
28#include <Singular/ipshell.h>
29#include <kernel/febase.h>
30#include <kernel/maps.h>
31#include <omalloc/omalloc.h>
32#include <kernel/kstd1.h>
33#include <kernel/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 (p!=NULL)
249        {
250          if( pIsConstant( p ) ) state= FglmHasOne;
251          else if ( (power= pIsPurePower( p )) > 0 )
252          {
253            fglmASSERT( 0 < power && power <= pVariables, "illegal power" );
254            if ( purePowers[power-1] == TRUE  ) state= FglmNotReduced;
255            else purePowers[power-1]= TRUE;
256          }
257          for ( int l = IDELEMS( theIdeal ) - 1; state == FglmOk && l >= 0; l-- )
258            if ( (k != l) && pDivisibleBy( p, (theIdeal->m)[l] ) )
259                state= FglmNotReduced;
260        }
261    }
262    if ( state == FglmOk )
263    {
264        for ( k= pVariables-1 ; (state == FglmOk) && (k >= 0); k-- )
265            if ( purePowers[k] == FALSE ) state= FglmNotZeroDim;
266    }
267    omFreeSize( (ADDRESS)purePowers, pVariables*sizeof( BOOLEAN ) );
268    return state;
269}
270
271// The main function for the fglm-Algorithm.
272// Checks the input-data, and calls fglmzero (see fglmzero.cc).
273// Returns the new groebnerbasis or 0 if an error occoured.
274BOOLEAN
275fglmProc( leftv result, leftv first, leftv second )
276{
277    FglmState state = FglmOk;
278
279    idhdl destRingHdl = currRingHdl;
280    ring destRing = currRing;
281    ideal destIdeal = NULL;
282    idhdl sourceRingHdl = (idhdl)first->data;
283    rSetHdl( sourceRingHdl );
284    ring sourceRing = currRing;
285
286    int * vperm = (int *)omAlloc0( (pVariables+1)*sizeof( int ) );
287    state= fglmConsistency( sourceRingHdl, destRingHdl, vperm );
288    omFreeSize( (ADDRESS)vperm, (pVariables+1)*sizeof(int) );
289
290    if ( state == FglmOk )
291    {
292        idhdl ih = currRing->idroot->get( second->Name(), myynest );
293        if ( (ih != NULL) && (IDTYP(ih)==IDEAL_CMD) )
294        {
295            ideal sourceIdeal;
296            if ( currQuotient != NULL )
297                sourceIdeal= fglmUpdatesource( IDIDEAL( ih ) );
298            else
299                sourceIdeal = IDIDEAL( ih );
300            state= fglmIdealcheck( sourceIdeal );
301            if ( state == FglmOk )
302            {
303                // Now the settings are compatible with FGLM
304                assumeStdFlag( (leftv)ih );
305                if ( fglmzero( IDRING(sourceRingHdl), sourceIdeal, destRingHdl, destIdeal, FALSE, (currQuotient != NULL) ) == FALSE )
306                    state= FglmNotReduced;
307            }
308        } else state= FglmNoIdeal;
309    }
310    if ( currRingHdl != destRingHdl )
311        rSetHdl( destRingHdl );
312    switch (state)
313    {
314        case FglmOk:
315            if ( currQuotient != NULL ) fglmUpdateresult( destIdeal );
316            break;
317        case FglmHasOne:
318            destIdeal= idInit(1,1);
319            (destIdeal->m)[0]= pOne();
320            state= FglmOk;
321            break;
322        case FglmIncompatibleRings:
323            Werror( "ring %s and current ring are incompatible", first->Name() );
324            destIdeal= idInit(0,0);
325            break;
326        case FglmNoIdeal:
327            Werror( "Can't find ideal %s in ring %s", second->Name(), first->Name() );
328            destIdeal= idInit(0,0);
329            break;
330        case FglmNotZeroDim:
331            Werror( "The ideal %s has to be 0-dimensional", second->Name() );
332            destIdeal= idInit(0,0);
333            break;
334        case FglmNotReduced:
335            Werror( "The ideal %s has to be given by a reduced SB", second->Name() );
336            destIdeal= idInit(0,0);
337            break;
338        default:
339            destIdeal= idInit(1,1);
340    }
341
342    result->rtyp = IDEAL_CMD;
343    result->data= (void *)destIdeal;
344    setFlag( result, FLAG_STD );
345    return (state != FglmOk);
346}
347
348// fglmQuotProc: Calculate I:f with FGLM methods.
349// Checks the input-data, and calls fglmquot (see fglmzero.cc).
350// Returns the new groebnerbasis if I:f or 0 if an error occoured.
351BOOLEAN
352fglmQuotProc( leftv result, leftv first, leftv second )
353{
354    FglmState state = FglmOk;
355
356    //    STICKYPROT("quotstart\n");
357    ideal sourceIdeal = (ideal)first->Data();
358    poly quot = (poly)second->Data();
359    ideal destIdeal = NULL;
360
361    state = fglmIdealcheck( sourceIdeal );
362    if ( state == FglmOk )
363    {
364      if ( quot == NULL ) state= FglmPolyIsZero;
365      else if ( pIsConstant( quot ) ) state= FglmPolyIsOne;
366    }
367
368    if ( state == FglmOk )
369    {
370      assumeStdFlag( first );
371      if ( fglmquot( sourceIdeal, quot, destIdeal ) == FALSE )
372        state= FglmNotReduced;
373    }
374
375    switch (state)
376    {
377        case FglmOk:
378            break;
379        case FglmHasOne:
380            destIdeal= idInit(1,1);
381            (destIdeal->m)[0]= pOne();
382            state= FglmOk;
383            break;
384        case FglmNotZeroDim:
385            Werror( "The ideal %s has to be 0-dimensional", first->Name() );
386            destIdeal= idInit(0,0);
387            break;
388        case FglmNotReduced:
389            Werror( "The poly %s has to be reduced", second->Name() );
390            destIdeal= idInit(0,0);
391            break;
392        case FglmPolyIsOne:
393            int k;
394            destIdeal= idInit( IDELEMS(sourceIdeal), 1 );
395            for ( k= IDELEMS( sourceIdeal )-1; k >=0; k-- )
396              (destIdeal->m)[k]= pCopy( (sourceIdeal->m)[k] );
397            state= FglmOk;
398            break;
399        case FglmPolyIsZero:
400            destIdeal= idInit(1,1);
401            (destIdeal->m)[0]= pOne();
402            state= FglmOk;
403            break;
404        default:
405            destIdeal= idInit(1,1);
406    }
407
408    result->rtyp = IDEAL_CMD;
409    result->data= (void *)destIdeal;
410    setFlag( result, FLAG_STD );
411    // STICKYPROT("quotend\n");
412    return (state != FglmOk);
413} // fglmQuotProt
414
415// The main function for finduni().
416// Checks the input-data, and calls FindUnivariateWrapper (see fglmzero.cc).
417// Returns an ideal containing the univariate Polynomials or 0 if an error
418// has occoured.
419BOOLEAN
420findUniProc( leftv result, leftv first )
421{
422    ideal sourceIdeal;
423    ideal destIdeal = NULL;
424    FglmState state;
425
426    sourceIdeal = (ideal)first->Data();
427
428    assumeStdFlag( first );
429    state= fglmIdealcheck( sourceIdeal );
430    if ( state == FglmOk )
431    {
432      // check for special cases: if the input contains
433      // univariate polys, try to reduce the problem
434      int i,k;
435      int count=0;
436      BOOLEAN * purePowers = (BOOLEAN *)omAlloc0( pVariables*sizeof( BOOLEAN ) );
437      for ( k= IDELEMS( sourceIdeal ) - 1; k >= 0; k-- )
438      {
439        if((i=pIsUnivariate(sourceIdeal->m[k]))>0)
440        {
441          if (purePowers[i-1]==0)
442          {
443            purePowers[i-1]=k;
444            count++;
445            if (count==pVariables) break;
446          }
447        }
448      }
449      if (count==pVariables)
450      {
451        destIdeal=idInit(pVariables,1);
452        for(k=pVariables-1; k>=0; k--) destIdeal->m[k]=pCopy(sourceIdeal->m[purePowers[k]]);
453      }
454      omFreeSize((ADDRESS)purePowers, pVariables*sizeof( BOOLEAN ) );
455      if (destIdeal!=NULL)
456            state = FglmOk;
457      else if ( FindUnivariateWrapper( sourceIdeal, destIdeal ) == FALSE )
458            state = FglmNotReduced;
459    }
460    switch (state)
461    {
462        case FglmOk:
463            break;
464        case FglmHasOne:
465            destIdeal= idInit(1,1);
466            (destIdeal->m)[0]= pOne();
467            state= FglmOk;
468            break;
469        case FglmNotZeroDim:
470            Werror( "The ideal %s has to be 0-dimensional", first->Name() );
471            destIdeal= idInit(0,0);
472            break;
473        case FglmNotReduced:
474            Werror( "The ideal %s has to be reduced", first->Name() );
475            destIdeal= idInit(0,0);
476            break;
477        default:
478            destIdeal= idInit(1,1);
479    }
480
481    result->rtyp = IDEAL_CMD;
482    result->data= (void *)destIdeal;
483
484    return FALSE;
485}
486#endif
487// ----------------------------------------------------------------------------
488// Local Variables: ***
489// compile-command: "make Singular" ***
490// page-delimiter: "^\\(\\|//!\\)" ***
491// fold-internal-margins: nil ***
492// End: ***
Note: See TracBrowser for help on using the repository browser.