source: git/kernel/fglmhom.cc @ 338842d

spielwiese
Last change on this file since 338842d was 35aab3, checked in by Hans Schönemann <hannes@…>, 20 years ago
This commit was generated by cvs2svn to compensate for changes in r6879, which included commits to RCS files with non-trunk default branches. git-svn-id: file:///usr/local/Singular/svn/trunk@6880 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.3 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fglmhom.cc,v 1.1.1.1 2003-10-06 12:15:53 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 "omalloc.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#ifndef HAVE_EXPLICIT_CONSTR
70    void initialize( poly m, int b, BOOLEAN ind )
71    {
72        basis = b;
73        inDest = ind;
74        mon.dm = m;
75        mon.sm = NULL;
76    }
77    void initialize( const homogElem h )
78    {
79        basis = h.basis;
80        inDest = h.inDest;
81        mon.dm = h.mon.dm;
82        mon.sm = h.mon.sm;
83    }
84#endif
85};
86
87struct homogData
88{
89    ideal sourceIdeal;
90    doublepoly * sourceHeads;
91    int numSourceHeads;
92    ideal destIdeal;
93    int numDestPolys;
94    homogElem * monlist;
95    int monlistmax;
96    int monlistblock;
97    int numMonoms;
98    int basisSize;
99    int overall;  // nur zum testen.
100    int numberofdestbasismonoms;
101//     homogData() : sourceHeads(NULL), numSourceHeads(0), monlist(NULL),
102//         numMonoms(0), basisSize(0) {}
103};
104
105int
106hfglmNextdegree( intvec * source, ideal current, int & deg )
107{
108    int numelems;
109    intvec * newhilb = hHstdSeries( current, NULL, currQuotient );
110
111    loop
112    {
113        if ( deg < newhilb->length() )
114        {
115            if ( deg < source->length() )
116                numelems= (*newhilb)[deg]-(*source)[deg];
117            else
118                numelems= (*newhilb)[deg];
119        }
120        else
121        {
122            if (deg < source->length())
123                numelems= -(*source)[deg];
124            else
125            {
126                deg= 0;
127                return 0;
128            }
129        }
130        if (numelems != 0)
131            return numelems;
132        deg++;
133    }
134    delete newhilb;
135}
136
137void
138generateMonoms( poly m, int var, int deg, homogData * dat )
139{
140    if ( var == pVariables ) {
141        BOOLEAN inSource = FALSE;
142        BOOLEAN inDest = FALSE;
143        poly mon = pCopy( m );
144        pSetExp( mon, var, deg );
145        pSetm( mon );
146        ++dat->overall;
147        int i;
148        for ( i= dat->numSourceHeads - 1; (i >= 0) && (inSource==FALSE); i-- ) {
149            if ( pDivisibleBy( dat->sourceHeads[i].dm, mon ) ) {
150                inSource= TRUE;
151            }
152        }
153        for ( i= dat->numDestPolys - 1; (i >= 0) && (inDest==FALSE); i-- ) {
154            if ( pDivisibleBy( (dat->destIdeal->m)[i], mon ) ) {
155                inDest= TRUE;
156            }
157        }
158        if ( (!inSource) || (!inDest) ) {
159            int basis = 0;
160            if ( !inSource )
161                basis= ++(dat->basisSize);
162            if ( !inDest )
163                ++dat->numberofdestbasismonoms;
164            if ( dat->numMonoms == dat->monlistmax ) {
165                int k;
166#ifdef HAVE_EXPLICIT_CONSTR
167                // Expand array using Singulars ReAlloc function
168                dat->monlist=
169                    (homogElem * )omReallocSize( dat->monlist,
170                                           (dat->monlistmax)*sizeof( homogElem ),
171                                           (dat->monlistmax+dat->monlistblock) * sizeof( homogElem ) );
172                for ( k= dat->monlistmax; k < (dat->monlistmax+dat->monlistblock); k++ )
173                    dat->monlist[k].homogElem();
174#else
175                // Expand array by generating new one and copying
176                int newsize = dat->monlistmax  + dat->monlistblock;
177                homogElem * tempelem = new homogElem[ newsize ];
178                // Copy old elements
179                for ( k= dat->monlistmax - 1; k >= 0; k-- )
180                    tempelem[k].initialize( dat->monlist[k] );
181                delete [] homogElem;
182                homogElem = tempelem;
183#endif
184                dat->monlistmax+= dat->monlistblock;
185            }
186#ifdef HAVE_EXPLICIT_CONSTR
187            dat->monlist[dat->numMonoms]= homogElem( mon, basis, inDest );
188#else
189            dat->monlist[dat->numMonoms].initialize( mon, basis, inDest );
190#endif
191            dat->numMonoms++;
192            if ( inSource && ! inDest ) PROT( "\\" );
193            if ( ! inSource && inDest ) PROT( "/" );
194            if ( ! inSource && ! inDest ) PROT( "." );
195        }
196        else {
197            pDelete( & mon );
198        }
199        return;
200    }
201    else {
202        poly newm = pCopy( m );
203        while ( deg >= 0 ) {
204            generateMonoms( newm, var+1, deg, dat );
205            pIncrExp( newm, var );
206            pSetm( newm );
207            deg--;
208        }
209        pDelete( & newm );
210    }
211    return;
212}
213
214void
215mapMonoms( ring oldRing, homogData & dat )
216{
217    int * vperm = (int *)omAlloc( (currRing->N + 1)*sizeof(int) );
218    maFindPerm( oldRing->names, oldRing->N, NULL, 0, currRing->names, currRing->N, NULL, 0, vperm, NULL );
219    //nSetMap( oldRing->ch, oldRing->parameter, oldRing->P, oldRing->minpoly );
220    nSetMap( oldRing );
221    int s;
222    for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
223//        dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, currRing->N, NULL, 0 );
224      // obachman: changed the folowing to reflect the new calling interface of
225      // pPermPoly -- Tim please check whether this is correct!
226        dat.monlist[s].mon.sm= pPermPoly( dat.monlist[s].mon.dm, vperm, oldRing, NULL, 0 );
227    }
228}
229
230void
231getVectorRep( homogData & dat )
232{
233    // Calculate the NormalForms
234    int s;
235    for ( s= 0;  s < dat.numMonoms; s++ ) {
236        if ( dat.monlist[s].inDest == FALSE ) {
237            fglmVector v;
238            if ( dat.monlist[s].basis == 0 ) {
239                v= fglmVector( dat.basisSize );
240                // now the monom is in L(source)
241                PROT( "(" );
242                poly nf = kNF( dat.sourceIdeal, NULL, dat.monlist[s].mon.sm );
243                PROT( ")" );
244                poly temp = nf;
245                while (temp != NULL ) {
246                    int t;
247                    for ( t= dat.numMonoms - 1; t >= 0; t-- ) {
248                        if ( dat.monlist[t].basis > 0 ) {
249                            if ( pLmEqual( dat.monlist[t].mon.sm, temp ) ) {
250                                number coeff= nCopy( pGetCoeff( temp ) );
251                                v.setelem( dat.monlist[t].basis, coeff );
252                            }
253                        }
254                    }
255                    pIter(temp);
256                }
257                pDelete( & nf );
258            }
259            else {
260                PROT( "." );
261                v= fglmVector( dat.basisSize, dat.monlist[s].basis );
262            }
263            dat.monlist[s].v= v;
264        }
265    }
266}
267
268void
269remapVectors( ring oldring, homogData & dat )
270{
271    //nSetMap( oldring->ch, oldring->parameter, oldring->P, oldring->minpoly );
272    nSetMap( oldring );
273    int s;
274    for ( s= dat.numMonoms - 1; s >= 0; s-- ) {
275        if ( dat.monlist[s].inDest == FALSE ) {
276            int k;
277            fglmVector newv( dat.basisSize );
278            for ( k= dat.basisSize; k > 0; k-- ){
279                number newnum= nMap( dat.monlist[s].v.getelem( k ) );
280                newv.setelem( k, newnum );
281            }
282            dat.monlist[s].dv= newv;
283        }
284    }
285}
286
287void
288gaussreduce( homogData & dat, int maxnum, int BS )
289{
290    int s;
291    int found= 0;
292
293    int destbasisSize = 0;
294    gaussReducer gauss( dat.basisSize );
295
296    for ( s= 0; (s < dat.numMonoms) && (found < maxnum); s++ ) {
297        if ( dat.monlist[s].inDest == FALSE ) {
298            if ( gauss.reduce( dat.monlist[s].dv ) == FALSE ) {
299                destbasisSize++;
300                dat.monlist[s].destbasis= destbasisSize;
301                gauss.store();
302                PROT( "." );
303            }
304            else {
305                fglmVector p= gauss.getDependence();
306                poly result = pCopy( dat.monlist[s].mon.dm );
307                pSetCoeff( result, nCopy( p.getconstelem( p.size() ) ) );
308                int l = 0;
309                int k;
310                for ( k= 1; k < p.size(); k++ ) {
311                    if ( ! p.elemIsZero( k ) ) {
312                        while ( dat.monlist[l].destbasis != k )
313                            l++;
314                        poly temp = pCopy( dat.monlist[l].mon.dm );
315                        pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
316                        result= pAdd( result, temp );
317                    }
318                }
319                if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
320//                PROT2( "(%s)", pString( result ) );
321                PROT( "+" );
322                found++;
323                (dat.destIdeal->m)[dat.numDestPolys]= result;
324                dat.numDestPolys++;
325                if ( IDELEMS(dat.destIdeal) == dat.numDestPolys ) {
326                    pEnlargeSet( & dat.destIdeal->m, IDELEMS( dat.destIdeal ), BS );
327                    IDELEMS( dat.destIdeal )+= BS;
328                }
329
330            }
331
332        }
333    }
334    PROT2( "(%i", s );
335    PROT2( "/%i)", dat.numberofdestbasismonoms );
336}
337
338
339BOOLEAN
340fglmhomog( idhdl sourceRingHdl, ideal sourceIdeal, idhdl destRingHdl, ideal & destIdeal )
341{
342#define groebnerBS 16
343    int numGBelems;
344    int deg = 0;
345
346    homogData dat;
347
348    // get the hilbert series and the leading monomials of the sourceIdeal:
349    rSetHdl( sourceRingHdl, TRUE );
350    ring sourceRing = currRing;
351
352    intvec * hilb = hHstdSeries( sourceIdeal, NULL, currQuotient );
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        dat.sourceHeads[s].sm= pHead( (sourceIdeal->m)[s] );
358    }
359    dat.numSourceHeads= IDELEMS( sourceIdeal );
360    rSetHdl( destRingHdl, TRUE );
361    ring destRing = currRing;
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        dat.sourceHeads[s].dm= pPermPoly( dat.sourceHeads[s].sm, vperm, sourceRing, NULL, 0 );
371    }
372
373    dat.destIdeal= idInit( groebnerBS, 1 );
374    dat.numDestPolys= 0;
375
376    while ( (numGBelems= hfglmNextdegree( hilb, dat.destIdeal, deg )) != 0 ) {
377        int num = 0;  // the number of monoms of degree deg
378        PROT2( "deg= %i ", deg );
379        PROT2( "num= %i\ngen>", numGBelems );
380        dat.monlistblock= 512;
381        dat.monlistmax= dat.monlistblock;
382#ifdef HAVE_EXPLICIT_CONSTR
383        dat.monlist= (homogElem *)omAlloc( dat.monlistmax*sizeof( homogElem ) );
384        int j;
385        for ( j= dat.monlistmax - 1; j >= 0; j-- ) dat.monlist[j].homogElem();
386#else
387        dat.monlist = new homogElem[ dat.monlistmax ];
388#endif
389        dat.numMonoms= 0;
390        dat.basisSize= 0;
391        dat.overall= 0;
392        dat.numberofdestbasismonoms= 0;
393
394        poly start= pOne();
395        generateMonoms( start, 1, deg, &dat );
396        pDelete( & start );
397
398        PROT2( "(%i/", dat.basisSize );
399        PROT2( "%i)\nvec>", dat.overall );
400        // switch to sourceRing and map monoms
401        rSetHdl( sourceRingHdl, TRUE );
402        mapMonoms( destRing, dat );
403        getVectorRep( dat );
404
405        // switch to destination Ring and remap the vectors
406        rSetHdl( destRingHdl, TRUE );
407        remapVectors( sourceRing, dat );
408
409        PROT( "<\nred>" );
410        // now do gaussian reduction
411        gaussreduce( dat, numGBelems, groebnerBS );
412
413#ifdef HAVE_EXPLICIT_CONSTR
414        omFreeSize( (ADDRESS)dat.monlist, dat.monlistmax*sizeof( homogElem ) );
415#else
416        delete [] dat.monlist;
417#endif
418        PROT( "<\n" );
419    }
420    PROT( "\n" );
421    destIdeal= dat.destIdeal;
422    idSkipZeroes( destIdeal );
423    return TRUE;
424}
425
426ideal
427fglmhomProc(leftv first, leftv second)
428{
429    idhdl dest= currRingHdl;
430    ideal result;
431    // in den durch das erste Argument angegeben Ring schalten:
432    rSetHdl( (idhdl)first->data, TRUE );
433    idhdl ih= currRing->idroot->get( second->name, myynest );
434    ASSERT( ih!=NULL, "Error: Can't find ideal in ring");
435    rSetHdl( dest, TRUE );
436
437    ideal i= IDIDEAL(ih);
438    fglmhomog( (idhdl)first->data, i, dest, result );
439
440    return( result );
441}
442
443#endif
444
445// Questions:
446// Muss ich einen intvec freigeben?
447// ----------------------------------------------------------------------------
448// Local Variables: ***
449// compile-command: "make Singular" ***
450// page-delimiter: "^\\(\\|//!\\)" ***
451// fold-internal-margins: nil ***
452// End: ***
Note: See TracBrowser for help on using the repository browser.