source: git/Singular/fglmcomb.cc @ 634dab0

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