source: git/kernel/fglm/fglmhom.cc @ 6cf934

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