source: git/Singular/fglmcomb.cc @ 48aa42

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