source: git/Singular/fglmhom.cc @ 5812c69

spielwiese
Last change on this file since 5812c69 was 5812c69, checked in by Hans Schönemann <hannes@…>, 26 years ago
cosmetic changes git-svn-id: file:///usr/local/Singular/svn/trunk@2517 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 12.2 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fglmhom.cc,v 1.11 1998-09-24 09:59:39 Singular Exp $
3
4/****************************************
5*  Computer Algebra System SINGULAR     *
6****************************************/
7/*
8* ABSTRACT - The FGLM-Algorithm extended for homogeneous ideals.
9*   Calculates via the hilbert-function a groebner basis.
10*/
11
12
13#include "mod2.h"
14#ifdef HAVE_FGLM_HOM
15#ifndef NOSTREAMIO
16#include <iostream.h>
17#endif
18#include "tok.h"
19#include "structs.h"
20#include "subexpr.h"
21#include "polys.h"
22#include "ideals.h"
23#include "ring.h"
24#include "ipid.h"
25#include "ipshell.h"
26#include "febase.h"
27#include "maps.h"
28#include "mmemory.h"
29#include "fglm.h"
30#include "fglmvec.h"
31#include "fglmgauss.h"
32#include "intvec.h"
33#include "kstd1.h"
34#include "stairc.h"  // -> hHStdSeries, hFirstSeries usw.
35#include <templates/ftmpl_list.h>
36
37// obachman: Got rid off those "redefiende messages by includeing fglm.h
38#include "fglm.h"
39#if 0
40#define PROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
41#define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
42#define PROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
43#define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
44#define fglmASSERT(ignore1,ignore2)
45#endif
46
47struct doublepoly
48{
49    poly sm;
50    poly dm;
51};
52
53class homogElem
54{
55public:
56    doublepoly mon;
57    fglmVector v;
58    fglmVector dv;
59    int basis;
60    int destbasis;
61    BOOLEAN inDest;
62    homogElem() : v(), dv(), basis(0), destbasis(0), inDest(FALSE) {}
63    homogElem( poly m, int b, BOOLEAN ind ) :
64        basis(b), inDest(ind)
65    {
66        mon.dm= m;
67        mon.sm= NULL;
68    }
69};
70
71struct homogData
72{
73    ideal sourceIdeal;
74    doublepoly * sourceHeads;
75    int numSourceHeads;
76    ideal destIdeal;
77    int numDestPolys;
78    homogElem * monlist;
79    int monlistmax;
80    int monlistblock;
81    int numMonoms;
82    int basisSize;
83    int overall;  // nur zum testen.
84    int numberofdestbasismonoms;
85//     homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL),
86//         numMonoms(0), basisSize(0) {}
87};
88
89int
90hfglmNextdegree( intvec * source, ideal current, int & deg )
91{
92    int numelems;
93    intvec * newhilb = hHstdSeries( current, NULL, currQuotient );
94
95    loop
96    {
97        if ( deg < newhilb->length() )
98        {
99            if ( deg < source->length() )
100                numelems= (*newhilb)[deg]-(*source)[deg];
101            else
102                numelems= (*newhilb)[deg];
103        }
104        else
105        {
106            if (deg < source->length())
107                numelems= -(*source)[deg];
108            else
109            {
110                deg= 0;
111                return 0;
112            }
113        }
114        if (numelems != 0)
115            return numelems;
116        deg++;
117    }
118    delete newhilb;
119}
120
121void
122generateMonoms( poly m, int var, int deg, homogData * dat )
123{
124    if ( var == pVariables ) {
125        BOOLEAN inSource = FALSE;
126        BOOLEAN inDest = FALSE;
127        poly mon = pCopy( m );
128        pSetExp( mon, var, deg );
129        pSetm( mon );
130        ++dat->overall;
131        int i;
132        for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) {
133            if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) {
134                inSource= TRUE;
135            }
136        }
137        for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) {
138            if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) {
139                inDest= TRUE;
140            }
141        }
142        if ( (!inSource) || (!inDest) ) {
143            int basis = 0;
144            if ( !inSource )
145                basis= ++(dat->basisSize);
146            if ( !inDest )
147                ++dat->numberofdestbasismonoms;
148            if ( dat->numMonoms == dat->monlistmax ) {
149                dat->monlist= (homogElem * )ReAlloc( dat->monlist, (dat->monlistmax)*sizeof( homogElem ), (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) );
150                int k;
151                for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
152                    dat->monlist[k].homogElem();
153                dat->monlistmax+= dat->monlistblock;
154            }
155            dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
156            dat->numMonoms++;
157            if ( inSource && ! inDest ) PROT( "\\" );
158            if ( ! inSource && inDest ) PROT( "/" );
159            if ( ! inSource && ! inDest ) PROT( "." );
160        }
161        else {
162            pDelete( & mon );
163        }
164        return;
165    }
166    else {
167        poly newm = pCopy( m );
168        while ( deg >= 0 ) {
169            generateMonoms( newm, var+1, deg, dat );
170            pIncrExp( newm, var );
171            pSetm( newm );
172            deg--;
173        }
174        pDelete( & newm );
175    }
176    return;
177}
178
179void
180mapMonoms( ring oldRing, homogData & dat )
181{
182    int * vperm = (int *)Alloc( (currRing->N + 1)*sizeof(int) );
183    maFindPerm( oldRing->names, oldRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL );
184    nSetMap( oldRing->ch, oldRing->parameter, oldRing->P, oldRing->minpoly );
185    int s;
186    for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
187//        dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 );
188      // obachman: changed the folowing to reflect the new calling interface of
189      // pPermPoly -- Tim please check whether this is correct!
190        dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 );
191    }
192}
193
194void
195getVectorRep( homogData & dat )
196{
197    // Calculate the NormalForms
198    int s;
199    for ( s= 0;  s < dat.numMonoms; s++ ) {
200        if ( dat.monlist[s].inDest == FALSE ) {
201            fglmVector v;
202            if ( dat.monlist[s].basis == 0 ) {
203                v= fglmVector( dat.basisSize );
204                // now the monom is in L(source)
205                PROT( "(" );
206                poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm );
207                PROT( ")" );
208                poly temp = nf;
209                while (temp != NULL ) {
210                    int t;
211                    for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
212                        if ( dat.monlist[t].basis > 0 ) {
213                            if ( pEqual( dat.monlist[t].mon.sm, temp ) ) {
214                                number coeff= nCopy( pGetCoeff( temp ) );
215                                v.setelem( dat.monlist[t].basis, coeff );
216                            }
217                        }
218                    }
219                    pIter(temp);
220                }
221                pDelete( & nf );
222            }
223            else {
224                PROT( "." );
225                v= fglmVector( dat.basisSize, dat.monlist[s].basis );
226            }
227            dat.monlist[s].v= v;
228        }
229    }
230}
231
232void
233remapVectors( ring oldring, homogData & dat )
234{
235    nSetMap( oldring->ch, oldring->parameter, oldring->P, oldring->minpoly );
236    int s;
237    for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
238        if ( dat.monlist[s].inDest == FALSE ) {
239            int k;
240            fglmVector newv( dat.basisSize );
241            for ( k= dat.basisSize; k > 0; k-- ){
242                number newnum= nMap( dat.monlist[s].v.getelem( k ) );
243                newv.setelem( k, newnum );
244            }
245            dat.monlist[s].dv= newv;
246        }
247    }
248}
249
250void
251gaussreduce( homogData & dat, int maxnum, int BS )
252{
253    int s;
254    int found= 0;
255
256    int destbasisSize = 0;
257    gaussReducer gauss( dat.basisSize );
258
259    for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) {
260        if ( dat.monlist[s].inDest == FALSE ) {
261            if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) {
262                destbasisSize++;
263                dat.monlist[s].destbasis= destbasisSize;
264                gauss.store();
265                PROT( "." );
266            }
267            else {
268                fglmVector p= gauss.getDependence();
269                poly result = pCopy( dat.monlist[s].mon.dm );
270                pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) );
271                int l = 0;
272                int k;
273                for ( k= 1; k < p.size(); k++ ) {
274                    if ( ! p.elemIsZero( k ) ) {
275                        while ( dat.monlist[l].destbasis != k )
276                            l++;
277                        poly temp = pCopy( dat.monlist[l].mon.dm );
278                        pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
279                        result= pAdd( result, temp );
280                    }
281                }
282                if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
283//                PROT2( "(%s)", pString( result ) );
284                PROT( "+" );
285                found++;
286                (dat.destIdeal->m)[dat.numDestPolys]= result;
287                dat.numDestPolys++;
288                if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
289                    pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS );
290                    IDELEMS( dat.destIdeal )+= BS;
291                }
292
293            }
294
295        }
296    }
297    PROT2( "(%i", s );
298    PROT2( "/%i)", dat.numberofdestbasismonoms );
299}
300
301
302BOOLEAN
303fglmhomog( idhdl sourceRingHdl, ideal sourceIdeal, idhdl destRingHdl, ideal & destIdeal )
304{
305#define groebnerBS 16
306    int numGBelems;
307    int deg = 0;
308
309    homogData dat;
310
311    // get the hilbert series and the leading monomials of the sourceIdeal:
312    rSetHdl( sourceRingHdl, TRUE );
313    ring sourceRing = currRing;
314
315    intvec * hilb = hHstdSeries( sourceIdeal, NULL, currQuotient );
316    int s;
317    dat.sourceIdeal= sourceIdeal;
318    dat.sourceHeads= (doublepoly *)Alloc( IDELEMS( sourceIdeal ) * sizeof( doublepoly ) );
319    for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- ) {
320        dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] );
321    }
322    dat.numSourceHeads= IDELEMS( sourceIdeal );
323    rSetHdl( destRingHdl, TRUE );
324    ring destRing = currRing;
325
326    // Map the sourceHeads to the destRing
327    int * vperm = (int *)Alloc( (sourceRing->N + 1)*sizeof(int) );
328    maFindPerm( sourceRing->names, sourceRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL );
329    nSetMap( sourceRing->ch, sourceRing->parameter, sourceRing->P, sourceRing->minpoly );
330    for ( s= IDELEMS( sourceIdeal ) - 1; s >= 0; s-- ) {
331        dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 );
332    }
333
334    dat.destIdeal= idInit( groebnerBS, 1 );
335    dat.numDestPolys= 0;
336
337    while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 ) {
338        int num = 0;  // the number of monoms of degree deg
339        PROT2( "deg= %i ", deg );
340        PROT2( "num= %i\ngen>", numGBelems );
341        dat.monlistblock= 512;
342        dat.monlistmax= dat.monlistblock;
343        dat.monlist= (homogElem *)Alloc( dat.monlistmax*sizeof( homogElem ) );
344        int j;
345        for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem();
346        dat.numMonoms= 0;
347        dat.basisSize= 0;
348        dat.overall= 0;
349        dat.numberofdestbasismonoms= 0;
350
351        poly start= pOne();
352        generateMonoms( start, 1, deg, &dat );
353        pDelete( & start );
354
355        PROT2( "(%i/", dat.basisSize );
356        PROT2( "%i)\nvec>", dat.overall );
357        // switch to sourceRing and map monoms
358        rSetHdl( sourceRingHdl, TRUE );
359        mapMonoms( destRing, dat );
360        getVectorRep( dat );
361
362        // switch to destination Ring and remap the vectors
363        rSetHdl( destRingHdl, TRUE );
364        remapVectors( sourceRing, dat );
365
366        PROT( "<\nred>" );
367        // now do gaussian reduction
368        gaussreduce( dat, numGBelems, groebnerBS );
369
370        Free( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) );
371        PROT( "<\n" );
372    }
373    PROT( "\n" );
374    destIdeal= dat.destIdeal;
375    idSkipZeroes( destIdeal );
376    return TRUE;
377}
378
379ideal
380fglmhomProc(leftv first, leftv second)
381{
382    idhdl dest= currRingHdl;
383    ideal result;
384    // in den durch das erste Argument angegeben Ring schalten:
385    rSetHdl( (idhdl)first->data, TRUE );
386    idhdl ih= currRing->idroot->get( second->name, myynest );
387    ASSERT( ih!=NULL, "Error: Can't find ideal in ring");
388    rSetHdl( dest, TRUE );
389
390    ideal i= IDIDEAL(ih);
391    fglmhomog( (idhdl)first->data, i, dest, result );
392
393    return( result );
394}
395
396#endif
397
398// Questions:
399// Muss ich einen intvec freigeben?
400// ----------------------------------------------------------------------------
401// Local Variables: ***
402// compile-command: "make Singular" ***
403// page-delimiter: "^\\(\\|//!\\)" ***
404// fold-internal-margins: nil ***
405// End: ***
Note: See TracBrowser for help on using the repository browser.