source: git/kernel/fglmcomb.cc @ fbc7cb

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