source: git/kernel/fglmhom.cc @ 762407

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