source: git/kernel/fglmcomb.cc @ abe5c8

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