source: git/kernel/fglmcomb.cc @ b29153

spielwiese
Last change on this file since b29153 was ea6e5e7, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: pSize ->polys.cc git-svn-id: file:///usr/local/Singular/svn/trunk@10243 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 15.4 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fglmcomb.cc,v 1.2 2007-07-24 15:10:51 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
34static void
35fglmEliminateMonomials( poly * pptr, fglmVector & v, polyset monomials, int numMonoms )
36{
37    poly temp = *pptr;
38    poly pretemp = NULL;
39    int point = 0;
40    int state;
41
42    while ( (temp != NULL) && (point < numMonoms) ) {
43        state= pCmp( temp, monomials[point] );
44        if ( state == 0 ) {
45            // Eliminate this monomial
46            poly todelete;
47            if ( pretemp == NULL ) {
48                todelete = temp;
49                pIter( *pptr );
50                temp= *pptr;
51            }
52            else {
53                todelete= temp;
54                pIter( temp );
55                pretemp->next= temp;
56            }
57            pGetCoeff( todelete )= nNeg( pGetCoeff( todelete ) );
58            number newelem = nAdd( pGetCoeff( todelete ), v.getconstelem( point+1 ) );
59            v.setelem( point+1, newelem );
60            nDelete( & pGetCoeff( todelete ) );
61            pLmFree( todelete );
62            point++;
63        }
64        else if ( state < 0 )
65            point++;
66        else {
67            pretemp= temp;
68            pIter( temp );
69        }
70    }
71}
72
73static BOOLEAN
74fglmReductionStep( poly * pptr, ideal source, int * w )
75{
76// returns TRUE if the leading monomial was reduced
77    if ( *pptr == NULL ) return FALSE;
78    int k;
79    int best = 0;
80    for ( k= IDELEMS( source ) - 1; k >= 0; k-- ) {
81        if ( pDivisibleBy( (source->m)[k], *pptr ) ) {
82            if ( best == 0 ) {
83                best= k + 1;
84            }
85            else {
86                if ( w[k] < w[best-1] ) {
87                    best= k + 1;
88                }
89            }
90        }
91    }
92    if ( best > 0 ) {
93        // OwnSPoly
94        poly p2 = (source->m)[best-1];
95        int i, diff;
96
97        poly m = pOne();
98        for ( i= pVariables; i > 0; i-- ) {
99            diff= pGetExp( *pptr, i ) - pGetExp( p2, i );
100            pSetExp( m, i, diff );
101        }
102        pSetm( m );
103        number n1 = nCopy( pGetCoeff( *pptr ) );
104        number n2 = pGetCoeff( p2 );
105
106        p2= pCopy( p2 );
107        pDeleteLm(pptr);
108        pDeleteLm( & p2 );
109        p2= pMult( m, p2 );
110
111        number temp = nDiv( n1, n2 );
112        n_Normalize( temp, currRing );
113        nDelete( & n1 );
114        n1= temp;
115        n1= nNeg( n1 );
116        pMult_nn( p2, n1 );
117//         pNormalize( p2 );
118        nDelete( & n1 );
119        *pptr= pAdd( *pptr, p2 );
120    }
121    return ( (best > 0) );
122}
123
124static void
125fglmReduce( poly * pptr, fglmVector & v, polyset m, int numMonoms, ideal source, int * w )
126{
127    BOOLEAN reduced = FALSE;
128    reduced= fglmReductionStep( pptr, source, w );
129    while ( reduced == TRUE ) {
130        fglmEliminateMonomials( pptr, v, m, numMonoms );
131        reduced= fglmReductionStep( pptr, source, w );
132    }
133    STICKYPROT( "<" );
134    poly temp = * pptr;
135    if ( temp != NULL ) {
136        while ( pNext( temp ) != NULL ) {
137            STICKYPROT( ">" );
138            reduced= fglmReductionStep( & pNext( temp ), source, w );
139            while ( reduced == TRUE ) {
140                fglmEliminateMonomials( & pNext( temp ), v, m, numMonoms );
141                reduced= fglmReductionStep( & pNext( temp ), source, w );
142            }
143            if ( pNext( temp ) != NULL ) {
144                pIter( temp );
145            }
146        }
147    }
148}
149
150static poly
151fglmNewLinearCombination( ideal source, poly monset )
152{
153    polyset m = NULL;
154    polyset nf = NULL;
155    fglmVector * mv = NULL;
156
157    fglmVector * v = NULL;
158    polyset basis = NULL;
159    int basisSize = 0;
160    int basisBS = 16;
161    int basisMax = basisBS;
162
163    int * weights = NULL;
164    int * lengthes = NULL;
165    int * order = NULL;
166
167    int numMonoms = 0;
168    int k;
169
170    // get number of monoms
171    numMonoms= pLength( monset );
172    STICKYPROT2( "%i monoms\n", numMonoms );
173
174    // Allcoate Memory and initialize sets
175    m= (polyset)omAlloc( numMonoms * sizeof( poly ) );
176    poly temp= monset;
177    for ( k= 0; k < numMonoms; k++ ) {
178//         m[k]= pOne();
179//         pSetExpV( m[k], temp->exp );
180//         pSetm( m[k] );
181        m[k]=pLmInit(temp);
182        pSetCoeff( m[k], nInit(1) );
183        pIter( temp );
184    }
185
186    nf= (polyset)omAlloc( numMonoms * sizeof( poly ) );
187
188#ifndef HAVE_EXPLICIT_CONSTR
189    mv= new fglmVector[ numMonoms ];
190#else
191    mv= (fglmVector *)omAlloc( numMonoms * sizeof( fglmVector ) );
192#endif
193
194#ifndef HAVE_EXPLICIT_CONSTR
195    v= new fglmVector[ numMonoms ];
196#else
197    v= (fglmVector *)omAlloc( numMonoms * sizeof( fglmVector ) );
198#endif
199
200    basisMax= basisBS;
201    basis= (polyset)omAlloc( basisMax * sizeof( poly ) );
202
203    weights= (int *)omAlloc( IDELEMS( source ) * sizeof( int ) );
204    STICKYPROT( "weights: " );
205    for ( k= 0; k < IDELEMS( source ); k++ ) {
206        poly temp= (source->m)[k];
207        int w = 0;
208        while ( temp != NULL ) {
209            w+= nSize( pGetCoeff( temp ) );
210            pIter( temp );
211        }
212        weights[k]= w;
213        STICKYPROT2( "%i ", w );
214    }
215    STICKYPROT( "\n" );
216    lengthes= (int *)omAlloc( numMonoms * sizeof( int ) );
217    order= (int *)omAlloc( numMonoms * sizeof( int ) );
218
219    // calculate the NormalForm in a special way
220    for ( k= 0; k < numMonoms; k++ )
221    {
222        STICKYPROT( "#" );
223        poly current = pCopy( m[k] );
224        fglmVector currV( numMonoms, k+1 );
225
226        fglmReduce( & current, currV, m, numMonoms, source, weights );
227        poly temp = current;
228        int b;
229        while ( temp != NULL )
230        {
231            BOOLEAN found = FALSE;
232            for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
233            {
234                if ( pLmEqual( temp, basis[b] ) )
235                {
236                    found= TRUE;
237                }
238            }
239            if ( found == FALSE )
240            {
241                if ( basisSize == basisMax )
242                {
243                    // Expand the basis
244                    basis= (polyset)omReallocSize( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
245                    basisMax+= basisBS;
246                }
247//                 basis[basisSize]= pOne();
248//                 pSetExpV( basis[basisSize], temp->exp );
249//                 pSetm( basis[basisSize] );
250                basis[basisSize]=pLmInit(temp);
251                pSetCoeff( basis[basisSize], nInit(1) );
252                basisSize++;
253            }
254            pIter( temp );
255        }
256        nf[k]= current;
257#ifndef HAVE_EXPLICIT_CONSTR
258        mv[k].mac_constr( currV );
259#else
260        mv[k].fglmVector( currV );
261#endif
262        STICKYPROT( "\n" );
263    }
264    // get the vector representation
265    for ( k= 0; k < numMonoms; k++ ) {
266        STICKYPROT( "." );
267
268#ifndef HAVE_EXPLICIT_CONSTR
269        v[k].mac_constr_i( basisSize );
270#else
271        v[k].fglmVector( basisSize );
272#endif
273        poly mon= nf[k];
274        while ( mon != NULL ) {
275            BOOLEAN found = FALSE;
276            int b= 0;
277            while ( found == FALSE ) {
278                if ( pLmEqual( mon, basis[b] ) ) {
279                    found= TRUE;
280                }
281                else {
282                    b++;
283                }
284            }
285            number coeff = nCopy( pGetCoeff( mon ) );
286            v[k].setelem( b+1, coeff );
287            pIter( mon );
288        }
289        pDelete( nf + k );
290    }
291    // Free Memory not needed anymore
292    omFreeSize( (ADDRESS)nf, numMonoms * sizeof( poly ) );
293    omFreeSize( (ADDRESS)weights, IDELEMS( source ) * sizeof( int ) );
294
295    STICKYPROT2( "\nbasis size: %i\n", basisSize );
296    STICKYPROT( "(clear basis" );
297    for ( k= 0; k < basisSize; k++ )
298        pDelete( basis + k );
299    STICKYPROT( ")\n" );
300    // gauss reduce
301    gaussReducer gauss( basisSize );
302    BOOLEAN isZero = FALSE;
303    fglmVector p;
304
305    STICKYPROT( "sizes: " );
306    for ( k= 0; k < numMonoms; k++ ) {
307        lengthes[k]= v[k].numNonZeroElems();
308        STICKYPROT2( "%i ", lengthes[k] );
309    }
310    STICKYPROT( "\n" );
311
312    int act = 0;
313    while ( (isZero == FALSE) && (act < numMonoms) ) {
314        int best = 0;
315        for ( k= numMonoms - 1; k >= 0; k-- ) {
316            if ( lengthes[k] > 0 ) {
317                if ( best == 0 ) {
318                    best= k+1;
319                }
320                else {
321                    if ( lengthes[k] < lengthes[best-1] ) {
322                        best= k+1;
323                    }
324                }
325            }
326        }
327        lengthes[best-1]= 0;
328        order[act]= best-1;
329        STICKYPROT2( " (%i) ", best );
330        if ( ( isZero= gauss.reduce( v[best-1] )) == TRUE ) {
331            p= gauss.getDependence();
332        }
333        else {
334            STICKYPROT( "+" );
335            gauss.store();
336            act++;
337        }
338#ifndef HAVE_EXPLICIT_CONSTR
339        v[best-1].clearelems();
340#else
341        v[best-1].~fglmVector();
342#endif
343    }
344    poly result = NULL;
345    if ( isZero == TRUE ) {
346        number gcd = p.gcd();
347        if ( ! nIsZero( gcd ) && ! ( nIsOne( gcd ) ) ) {
348            p/= gcd;
349        }
350        nDelete( & gcd );
351        fglmVector temp( numMonoms );
352        for ( k= 0; k < p.size(); k++ ) {
353            if ( ! p.elemIsZero( k+1 ) ) {
354                temp+= p.getconstelem( k+1 ) * mv[order[k]];
355            }
356        }
357        number denom = temp.clearDenom();
358        nDelete( & denom );
359        gcd = temp.gcd();
360        if ( ! nIsZero( gcd ) && ! nIsOne( gcd ) ) {
361            temp/= gcd;
362        }
363        nDelete( & gcd );
364
365        poly sum = NULL;
366        for ( k= 1; k <= numMonoms; k++ ) {
367            if ( ! temp.elemIsZero( k ) ) {
368                if ( result == NULL ) {
369                    result= pCopy( m[k-1] );
370                    sum= result;
371                }
372                else {
373                    sum->next= pCopy( m[k-1] );
374                    pIter( sum );
375                }
376                pSetCoeff( sum, nCopy( temp.getconstelem( k ) ) );
377            }
378        }
379        pContent( result );
380        if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
381    }
382    // Free Memory
383    omFreeSize( (ADDRESS)lengthes, numMonoms * sizeof( int ) );
384    omFreeSize( (ADDRESS)order, numMonoms * sizeof( int ) );
385//     for ( k= 0; k < numMonoms; k++ )
386//         v[k].~fglmVector();
387#ifndef HAVE_EXPLICIT_CONSTR
388    delete [] v;
389#else
390    omFreeSize( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
391#endif
392
393    for ( k= 0; k < basisSize; k++ )
394        pDelete( basis + k );
395    omFreeSize( (ADDRESS)basis, basisMax * sizeof( poly ) );
396
397#ifndef HAVE_EXPLICIT_CONSTR
398    delete [] mv;
399#else
400    for ( k= 0; k < numMonoms; k++ )
401        mv[k].~fglmVector();
402    omFreeSize( (ADDRESS)mv, numMonoms * sizeof( fglmVector ) );
403#endif
404
405    for ( k= 0; k < numMonoms; k++ )
406        pDelete( m + k );
407    omFreeSize( (ADDRESS)m, numMonoms * sizeof( poly ) );
408
409    STICKYPROT( "\n" );
410    return result;
411}
412
413
414static poly
415fglmLinearCombination( 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#endif
551// Local Variables: ***
552// compile-command: "make Singular" ***
553// page-delimiter: "^\\(\\|//!\\)" ***
554// fold-internal-margins: nil ***
555// End: ***
Note: See TracBrowser for help on using the repository browser.