source: git/Singular/fglmcomb.cc @ 3bf67ab

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