source: git/kernel/fglmhom.cc @ fbc7cb

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