source: git/kernel/fglmcomb.cc @ 68349d

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