source: git/Singular/fglmcomb.cc @ cf80ce

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