source: git/kernel/fglm/fglmcomb.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: 15.3 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2
3/****************************************
4*  Computer Algebra System SINGULAR     *
5****************************************/
6/*
7* ABSTRACT -
8*/
9
10
11
12
13#include "kernel/mod2.h"
14
15#include "factory/factory.h"
16#include "misc/options.h"
17#include "kernel/polys.h"
18#include "kernel/ideals.h"
19#include "polys/monomials/ring.h"
20#include "polys/monomials/maps.h"
21#include "fglmvec.h"
22#include "fglmgauss.h"
23#include "kernel/GBEngine/kstd1.h"
24
25#include "fglm.h"
26
27#ifndef NOSTREAMIO
28#  ifdef HAVE_IOSTREAM
29#    include <iostream>
30#  else
31#    include <iostream.h>
32#  endif
33#endif
34
35static void
36fglmEliminateMonomials( poly * pptr, fglmVector & v, polyset monomials, int numMonoms )
37{
38    poly temp = *pptr;
39    poly pretemp = NULL;
40    int point = 0;
41    int state;
42
43    while ( (temp != NULL) && (point < numMonoms) ) {
44        state= pCmp( temp, monomials[point] );
45        if ( state == 0 ) {
46            // Eliminate this monomial
47            poly todelete;
48            if ( pretemp == NULL ) {
49                todelete = temp;
50                pIter( *pptr );
51                temp= *pptr;
52            }
53            else {
54                todelete= temp;
55                pIter( temp );
56                pretemp->next= temp;
57            }
58            pGetCoeff( todelete )= nInpNeg( pGetCoeff( todelete ) );
59            number newelem = nAdd( pGetCoeff( todelete ), v.getconstelem( point+1 ) );
60            v.setelem( point+1, newelem );
61            nDelete( & pGetCoeff( todelete ) );
62            pLmFree( todelete );
63            point++;
64        }
65        else if ( state < 0 )
66            point++;
67        else {
68            pretemp= temp;
69            pIter( temp );
70        }
71    }
72}
73
74static BOOLEAN
75fglmReductionStep( poly * pptr, ideal source, int * w )
76{
77// returns TRUE if the leading monomial was reduced
78    if ( *pptr == NULL ) return FALSE;
79    int k;
80    int best = 0;
81    for ( k= IDELEMS( source ) - 1; k >= 0; k-- ) {
82        if ( pDivisibleBy( (source->m)[k], *pptr ) ) {
83            if ( best == 0 ) {
84                best= k + 1;
85            }
86            else {
87                if ( w[k] < w[best-1] ) {
88                    best= k + 1;
89                }
90            }
91        }
92    }
93    if ( best > 0 )
94    {
95        // OwnSPoly
96        poly p2 = (source->m)[best-1];
97        int i, diff;
98
99        poly m = pOne();
100        for ( i= (currRing->N); i > 0; i-- )
101        {
102            diff= pGetExp( *pptr, i ) - pGetExp( p2, i );
103            pSetExp( m, i, diff );
104        }
105        pSetm( m );
106        number n1 = nCopy( pGetCoeff( *pptr ) );
107        number n2 = pGetCoeff( p2 );
108
109        p2= pCopy( p2 );
110        pLmDelete(pptr);
111        pLmDelete( & p2 );
112        p2= pMult( m, p2 );
113
114        number temp = nDiv( n1, n2 );
115        n_Normalize( temp, currRing->cf );
116        nDelete( & n1 );
117        n1= temp;
118        n1= nInpNeg( n1 );
119        p2=__p_Mult_nn( p2, n1, currRing );
120//         pNormalize( p2 );
121        nDelete( & n1 );
122        *pptr= pAdd( *pptr, p2 );
123    }
124    return ( (best > 0) );
125}
126
127static void
128fglmReduce( poly * pptr, fglmVector & v, polyset m, int numMonoms, ideal source, int * w )
129{
130    BOOLEAN reduced = FALSE;
131    reduced= fglmReductionStep( pptr, source, w );
132    while ( reduced == TRUE ) {
133        fglmEliminateMonomials( pptr, v, m, numMonoms );
134        reduced= fglmReductionStep( pptr, source, w );
135    }
136    STICKYPROT( "<" );
137    poly temp = * pptr;
138    if ( temp != NULL ) {
139        while ( pNext( temp ) != NULL ) {
140            STICKYPROT( ">" );
141            reduced= fglmReductionStep( & pNext( temp ), source, w );
142            while ( reduced == TRUE ) {
143                fglmEliminateMonomials( & pNext( temp ), v, m, numMonoms );
144                reduced= fglmReductionStep( & pNext( temp ), source, w );
145            }
146            if ( pNext( temp ) != NULL ) {
147                pIter( temp );
148            }
149        }
150    }
151}
152
153poly fglmNewLinearCombination( ideal source, poly monset )
154{
155    polyset m = NULL;
156    polyset nf = NULL;
157    fglmVector * mv = NULL;
158
159    fglmVector * v = NULL;
160    polyset basis = NULL;
161    int basisSize = 0;
162    int basisBS = 16;
163    int basisMax = basisBS;
164
165    int * weights = NULL;
166    int * lengths = NULL;
167    int * order = NULL;
168
169    int numMonoms = 0;
170    int k;
171
172    // get number of monoms
173    numMonoms= pLength( monset );
174    STICKYPROT2( "%i monoms\n", numMonoms );
175
176    // Allocate Memory and initialize sets
177    m= (polyset)omAlloc( numMonoms * sizeof( poly ) );
178    poly temp= monset;
179    for ( k= 0; k < numMonoms; k++ ) {
180//         m[k]= pOne();
181//         pSetExpV( m[k], temp->exp );
182//         pSetm( m[k] );
183        m[k]=pLmInit(temp);
184        pSetCoeff( m[k], nInit(1) );
185        pIter( temp );
186    }
187
188    nf= (polyset)omAlloc( numMonoms * sizeof( poly ) );
189
190#ifndef HAVE_EXPLICIT_CONSTR
191    mv= new fglmVector[ numMonoms ];
192#else
193    mv= (fglmVector *)omAlloc( numMonoms * sizeof( fglmVector ) );
194#endif
195
196#ifndef HAVE_EXPLICIT_CONSTR
197    v= new fglmVector[ numMonoms ];
198#else
199    v= (fglmVector *)omAlloc( numMonoms * sizeof( fglmVector ) );
200#endif
201
202    basisMax= basisBS;
203    basis= (polyset)omAlloc( basisMax * sizeof( poly ) );
204
205    weights= (int *)omAlloc( IDELEMS( source ) * sizeof( int ) );
206    STICKYPROT( "weights: " );
207    for ( k= 0; k < IDELEMS( source ); k++ ) {
208        poly temp= (source->m)[k];
209        int w = 0;
210        while ( temp != NULL ) {
211            w+= nSize( pGetCoeff( temp ) );
212            pIter( temp );
213        }
214        weights[k]= w;
215        STICKYPROT2( "%i ", w );
216    }
217    STICKYPROT( "\n" );
218    lengths= (int *)omAlloc( numMonoms * sizeof( int ) );
219    order= (int *)omAlloc( numMonoms * sizeof( int ) );
220
221    // calculate the NormalForm in a special way
222    for ( k= 0; k < numMonoms; k++ )
223    {
224        STICKYPROT( "#" );
225        poly current = pCopy( m[k] );
226        fglmVector currV( numMonoms, k+1 );
227
228        fglmReduce( & current, currV, m, numMonoms, source, weights );
229        poly temp = current;
230        int b;
231        while ( temp != NULL )
232        {
233            BOOLEAN found = FALSE;
234            for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
235            {
236                if ( pLmEqual( temp, basis[b] ) )
237                {
238                    found= TRUE;
239                }
240            }
241            if ( found == FALSE )
242            {
243                if ( basisSize == basisMax )
244                {
245                    // Expand the basis
246                    basis= (polyset)omReallocSize( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
247                    basisMax+= basisBS;
248                }
249//                 basis[basisSize]= pOne();
250//                 pSetExpV( basis[basisSize], temp->exp );
251//                 pSetm( basis[basisSize] );
252                basis[basisSize]=pLmInit(temp);
253                pSetCoeff( basis[basisSize], nInit(1) );
254                basisSize++;
255            }
256            pIter( temp );
257        }
258        nf[k]= current;
259#ifndef HAVE_EXPLICIT_CONSTR
260        mv[k].mac_constr( currV );
261#else
262        mv[k].fglmVector( currV );
263#endif
264        STICKYPROT( "\n" );
265    }
266    // get the vector representation
267    for ( k= 0; k < numMonoms; k++ ) {
268        STICKYPROT( "." );
269
270#ifndef HAVE_EXPLICIT_CONSTR
271        v[k].mac_constr_i( basisSize );
272#else
273        v[k].fglmVector( basisSize );
274#endif
275        poly mon= nf[k];
276        while ( mon != NULL ) {
277            BOOLEAN found = FALSE;
278            int b= 0;
279            while ( found == FALSE ) {
280                if ( pLmEqual( mon, basis[b] ) ) {
281                    found= TRUE;
282                }
283                else {
284                    b++;
285                }
286            }
287            number coeff = nCopy( pGetCoeff( mon ) );
288            v[k].setelem( b+1, coeff );
289            pIter( mon );
290        }
291        pDelete( nf + k );
292    }
293    // Free Memory not needed anymore
294    omFreeSize( (ADDRESS)nf, numMonoms * sizeof( poly ) );
295    omFreeSize( (ADDRESS)weights, IDELEMS( source ) * sizeof( int ) );
296
297    STICKYPROT2( "\nbasis size: %i\n", basisSize );
298    STICKYPROT( "(clear basis" );
299    for ( k= 0; k < basisSize; k++ )
300        pDelete( basis + k );
301    STICKYPROT( ")\n" );
302    // gauss reduce
303    gaussReducer gauss( basisSize );
304    BOOLEAN isZero = FALSE;
305    fglmVector p;
306
307    STICKYPROT( "sizes: " );
308    for ( k= 0; k < numMonoms; k++ ) {
309        lengths[k]= v[k].numNonZeroElems();
310        STICKYPROT2( "%i ", lengths[k] );
311    }
312    STICKYPROT( "\n" );
313
314    int act = 0;
315    while ( (isZero == FALSE) && (act < numMonoms) ) {
316        int best = 0;
317        for ( k= numMonoms - 1; k >= 0; k-- ) {
318            if ( lengths[k] > 0 ) {
319                if ( best == 0 ) {
320                    best= k+1;
321                }
322                else {
323                    if ( lengths[k] < lengths[best-1] ) {
324                        best= k+1;
325                    }
326                }
327            }
328        }
329        lengths[best-1]= 0;
330        order[act]= best-1;
331        STICKYPROT2( " (%i) ", best );
332        if ( ( isZero= gauss.reduce( v[best-1] )) == TRUE ) {
333            p= gauss.getDependence();
334        }
335        else {
336            STICKYPROT( "+" );
337            gauss.store();
338            act++;
339        }
340#ifndef HAVE_EXPLICIT_CONSTR
341        v[best-1].clearelems();
342#else
343        v[best-1].~fglmVector();
344#endif
345    }
346    poly result = NULL;
347    if ( isZero == TRUE ) {
348        number gcd = p.gcd();
349        if ( ! nIsZero( gcd ) && ! ( nIsOne( gcd ) ) ) {
350            p/= gcd;
351        }
352        nDelete( & gcd );
353        fglmVector temp( numMonoms );
354        for ( k= 0; k < p.size(); k++ ) {
355            if ( ! p.elemIsZero( k+1 ) ) {
356                temp+= p.getconstelem( k+1 ) * mv[order[k]];
357            }
358        }
359        number denom = temp.clearDenom();
360        nDelete( & denom );
361        gcd = temp.gcd();
362        if ( ! nIsZero( gcd ) && ! nIsOne( gcd ) ) {
363            temp/= gcd;
364        }
365        nDelete( & gcd );
366
367        poly sum = NULL;
368        for ( k= 1; k <= numMonoms; k++ ) {
369            if ( ! temp.elemIsZero( k ) ) {
370                if ( result == NULL ) {
371                    result= pCopy( m[k-1] );
372                    sum= result;
373                }
374                else {
375                    sum->next= pCopy( m[k-1] );
376                    pIter( sum );
377                }
378                pSetCoeff( sum, nCopy( temp.getconstelem( k ) ) );
379            }
380        }
381        p_Cleardenom( result, currRing );
382    }
383    // Free Memory
384    omFreeSize( (ADDRESS)lengths, numMonoms * sizeof( int ) );
385    omFreeSize( (ADDRESS)order, numMonoms * sizeof( int ) );
386//     for ( k= 0; k < numMonoms; k++ )
387//         v[k].~fglmVector();
388#ifndef HAVE_EXPLICIT_CONSTR
389    delete [] v;
390#else
391    omFreeSize( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
392#endif
393
394    for ( k= 0; k < basisSize; k++ )
395        pDelete( basis + k );
396    omFreeSize( (ADDRESS)basis, basisMax * sizeof( poly ) );
397
398#ifndef HAVE_EXPLICIT_CONSTR
399    delete [] mv;
400#else
401    for ( k= 0; k < numMonoms; k++ )
402        mv[k].~fglmVector();
403    omFreeSize( (ADDRESS)mv, numMonoms * sizeof( fglmVector ) );
404#endif
405
406    for ( k= 0; k < numMonoms; k++ )
407        pDelete( m + k );
408    omFreeSize( (ADDRESS)m, numMonoms * sizeof( poly ) );
409
410    STICKYPROT( "\n" );
411    return result;
412}
413
414
415poly fglmLinearCombination( ideal source, poly monset )
416{
417    int k;
418    poly temp;
419
420    polyset m;
421    polyset nf;
422    fglmVector * v;
423
424    polyset basis;
425    int basisSize = 0;
426    int basisBS = 16;
427    int basisMax;
428    // get number of monomials in monset
429    int numMonoms = 0;
430    temp = monset;
431    while ( temp != NULL ) {
432        numMonoms++;
433        pIter( temp );
434    }
435    // Allocate Memory
436    m= (polyset)omAlloc( numMonoms * sizeof( poly ) );
437    nf= (polyset)omAlloc( numMonoms * sizeof( poly ) );
438    // Warning: The fglmVectors in v are yet not initialized
439    v= (fglmVector *)omAlloc( numMonoms * sizeof( fglmVector ) );
440    basisMax= basisBS;
441    basis= (polyset)omAlloc( basisMax * sizeof( poly ) );
442
443    // get the NormalForm and the basis monomials
444    temp= monset;
445    for ( k= 0; k < numMonoms; k++ ) {
446        poly mon= pHead( temp );
447        if ( ! nIsOne( pGetCoeff( mon ) ) ) {
448            nDelete( & pGetCoeff( mon ) );
449            pSetCoeff( mon, nInit( 1 ) );
450        }
451        STICKYPROT( "(" );
452        nf[k]= kNF( source, currRing->qideal, mon );
453        STICKYPROT( ")" );
454
455        // search through basis
456        STICKYPROT( "[" );
457        poly sm = nf[k];
458        while ( sm != NULL ) {
459            BOOLEAN found = FALSE;
460            int b;
461            for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
462                if ( pLmEqual( sm, basis[b] ) ) found= TRUE;
463            if ( found == FALSE ) {
464                // Expand the basis
465                if ( basisSize == basisMax ) {
466                    basis= (polyset)omReallocSize( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
467                    basisMax+= basisBS;
468                }
469                basis[basisSize]= pHead( sm );
470                nDelete( & pGetCoeff( basis[basisSize] ) );
471                pSetCoeff( basis[basisSize], nInit( 1 ) );
472                basisSize++;
473            }
474            pIter( sm );
475        }
476        STICKYPROT( "]" );
477        m[k]= mon;
478        pIter( temp );
479    }
480    // get the vector representation
481    STICKYPROT2( "(%i)", basisSize );
482    for ( k= 0; k < numMonoms; k++ ) {
483#ifndef HAVE_EXPLICIT_CONSTR
484        v[k].mac_constr_i( basisSize );
485#else
486        v[k].fglmVector( basisSize );
487#endif
488        STICKYPROT( "(+" );
489        poly mon= nf[k];
490        while ( mon != NULL ) {
491            BOOLEAN found = FALSE;
492            int b= 0;
493            while ( found == FALSE ) {
494                if ( pLmEqual( mon, basis[b] ) )
495                    found= TRUE;
496                else
497                    b++;
498            }
499            number coeff = nCopy( pGetCoeff( mon ) );
500            v[k].setelem( b+1, coeff );
501            pIter( mon );
502        }
503        STICKYPROT( ")" );
504    }
505    // gauss reduce
506    gaussReducer gauss( basisSize );
507    BOOLEAN isZero = FALSE;
508    fglmVector p;
509    for ( k= 0; (k < numMonoms) && (isZero == FALSE); k++ ) {
510        STICKYPROT( "(-" );
511        if ( ( isZero= gauss.reduce( v[k] )) == TRUE )
512            p= gauss.getDependence();
513        else
514            gauss.store();
515        STICKYPROT( ")" );
516    }
517    poly comb = NULL;
518    if ( isZero == TRUE ) {
519        number gcd = p.gcd();
520        if ( ! nIsZero( gcd ) && ! ( nIsOne( gcd ) ) )
521            p/= gcd;
522        nDelete( & gcd );
523        for ( k= 1; k <= p.size(); k++ ) {
524            if ( ! p.elemIsZero( k ) ) {
525                poly temp = pCopy( m[k-1] );
526                pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
527                comb= pAdd( comb, temp );
528            }
529        }
530        if ( ! nGreaterZero( pGetCoeff( comb ) ) ) comb= pNeg( comb );
531    }
532
533    // Free Memory
534    for ( k= 0; k < numMonoms; k++ ) {
535        pDelete( m + k );
536        pDelete( nf + k );
537    }
538    omFreeSize( (ADDRESS)m, numMonoms * sizeof( poly ) );
539    omFreeSize( (ADDRESS)nf, numMonoms * sizeof( poly ) );
540    // Warning: At this point all Vectors in v have to be initialized
541    for ( k= numMonoms - 1; k >= 0; k-- ) v[k].~fglmVector();
542    omFreeSize( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
543    for ( k= 0; k < basisSize; k++ )
544        pDelete( basis + k );
545    omFreeSize( (ADDRESS)basis, basisMax * sizeof( poly ) );
546    STICKYPROT( "\n" );
547    return comb;
548}
549
550// Local Variables: ***
551// compile-command: "make Singular" ***
552// page-delimiter: "^\\(\\|//!\\)" ***
553// fold-internal-margins: nil ***
554// End: ***
Note: See TracBrowser for help on using the repository browser.