source: git/kernel/fglmcomb.cc @ 599326

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