source: git/Singular/fglmcomb.cc @ c4bbf1f

spielwiese
Last change on this file since c4bbf1f was e3cab3, checked in by Wilfred Pohl <pohl@…>, 26 years ago
substitute macintosh by __MWERKS__ git-svn-id: file:///usr/local/Singular/svn/trunk@1359 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.9 1998-04-08 12:11:24 pohl 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/ftmpl_list.h>
34
35// nur fuer debug-Ausgaben:
36static int
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
47static void
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
86static BOOLEAN
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
137static void
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
163static poly
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        m[k]= pNew();
195        pCopy2( m[k], temp );
196        pSetCoeff( m[k], nInit(1) );
197        temp= pIter( temp );
198    }
199
200    nf= (polyset)Alloc( numMonoms * sizeof( poly ) );
201    mv= (fglmVector *)Alloc( numMonoms * sizeof( fglmVector ) );
202
203    v= (fglmVector *)Alloc( numMonoms * sizeof( fglmVector ) );
204    basisMax= basisBS;
205    basis= (polyset)Alloc( basisMax * sizeof( poly ) );
206
207    weights= (int *)Alloc( IDELEMS( source ) * sizeof( int ) );
208    STICKYPROT( "weights: " );
209    for ( k= 0; k < IDELEMS( source ); k++ ) {
210        poly temp= (source->m)[k];
211        int w = 0;
212        while ( temp != NULL ) {
213            w+= nSize( pGetCoeff( temp ) );
214            temp= pIter( temp );
215        }
216        weights[k]= w;
217        STICKYPROT2( "%i ", w );
218    }
219    STICKYPROT( "\n" );
220    lengthes= (int *)Alloc( numMonoms * sizeof( int ) );
221    order= (int *)Alloc( numMonoms * sizeof( int ) );
222
223    // calculate the NormalForm in a special way
224    for ( k= 0; k < numMonoms; k++ )
225    {
226        STICKYPROT( "#" );
227        poly current = pCopy( m[k] );
228        fglmVector currV( numMonoms, k+1 );
229
230        fglmReduce( & current, currV, m, numMonoms, source, weights );
231        poly temp = current;
232        int b;
233        while ( temp != NULL )
234        {
235            BOOLEAN found = FALSE;
236            for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
237            {
238                if ( pEqual( temp, basis[b] ) )
239                {
240                    found= TRUE;
241                }
242            }
243            if ( found == FALSE )
244            {
245                if ( basisSize == basisMax )
246                {
247                    // Expand the basis
248                    basis= (polyset)ReAlloc( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
249                    basisMax+= basisBS;
250                }
251//                 basis[basisSize]= pOne();
252//                 pSetExpV( basis[basisSize], temp->exp );
253//                 pSetm( basis[basisSize] );
254                basis[basisSize]= pNew();
255                pCopy2( basis[basisSize], temp );
256                pSetCoeff( basis[basisSize], nInit(1) );
257                basisSize++;
258            }
259            temp= pIter( temp );
260        }
261        nf[k]= current;
262#ifdef __MWERKS__
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#ifdef __MWERKS__
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 ( pEqual( 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            mon= pIter( mon );
293        }
294        pDelete( nf + k );
295    }
296    // Free Memory not needed anymore
297    Free( (ADDRESS)nf, numMonoms * sizeof( poly ) );
298    Free( (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        v[best-1].~fglmVector();
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        pContent( result );
381        if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
382    }
383    // Free Memory
384    Free( (ADDRESS)lengthes, numMonoms * sizeof( int ) );
385    Free( (ADDRESS)order, numMonoms * sizeof( int ) );
386//     for ( k= 0; k < numMonoms; k++ )
387//         v[k].~fglmVector();
388    Free( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
389    for ( k= 0; k < basisSize; k++ )
390        pDelete( basis + k );
391    Free( (ADDRESS)basis, basisMax * sizeof( poly ) );
392    for ( k= 0; k < numMonoms; k++ )
393        mv[k].~fglmVector();
394    Free( (ADDRESS)mv, numMonoms * sizeof( fglmVector ) );
395
396    for ( k= 0; k < numMonoms; k++ )
397        pDelete( m + k );
398    Free( (ADDRESS)m, numMonoms * sizeof( poly ) );
399
400    STICKYPROT( "\n" );
401    return result;
402}
403
404
405static poly
406fglmLinearCombination( ideal source, poly monset )
407{
408    int k;
409    poly temp;
410
411    polyset m;
412    polyset nf;
413    fglmVector * v;
414
415    polyset basis;
416    int basisSize = 0;
417    int basisBS = 16;
418    int basisMax;
419    // get number of monomials in monset
420    int numMonoms = 0;
421    temp = monset;
422    while ( temp != NULL ) {
423        numMonoms++;
424        temp= pIter( temp );
425    }
426    // Allocate Memory
427    m= (polyset)Alloc( numMonoms * sizeof( poly ) );
428    nf= (polyset)Alloc( numMonoms * sizeof( poly ) );
429    // Warning: The fglmVectors in v are yet not initialized
430    v= (fglmVector *)Alloc( numMonoms * sizeof( fglmVector ) );
431    basisMax= basisBS;
432    basis= (polyset)Alloc( basisMax * sizeof( poly ) );
433
434    // get the NormalForm and the basis monomials
435    temp= monset;
436    for ( k= 0; k < numMonoms; k++ ) {
437        poly mon= pHead( temp );
438        if ( ! nIsOne( pGetCoeff( mon ) ) ) {
439            nDelete( & pGetCoeff( mon ) );
440            pSetCoeff( mon, nInit( 1 ) );
441        }
442        STICKYPROT( "(" );
443        nf[k]= kNF( source, currRing->qideal, mon );
444        STICKYPROT( ")" );
445
446        // search through basis
447        STICKYPROT( "[" );
448        poly sm = nf[k];
449        while ( sm != NULL ) {
450            BOOLEAN found = FALSE;
451            int b;
452            for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
453                if ( pEqual( sm, basis[b] ) ) found= TRUE;
454            if ( found == FALSE ) {
455                // Expand the basis
456                if ( basisSize == basisMax ) {
457                    basis= (polyset)ReAlloc( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
458                    basisMax+= basisBS;
459                }
460                basis[basisSize]= pHead( sm );
461                nDelete( & pGetCoeff( basis[basisSize] ) );
462                pSetCoeff( basis[basisSize], nInit( 1 ) );
463                basisSize++;
464            }
465            sm= pIter( sm );
466        }
467        STICKYPROT( "]" );
468        m[k]= mon;
469        temp= pIter( temp );
470    }
471    // get the vector representation
472    STICKYPROT2( "(%i)", basisSize );
473    for ( k= 0; k < numMonoms; k++ ) {
474#ifdef __MWERKS__
475        v[k].mac_constr_i( basisSize );
476#else
477        v[k].fglmVector( basisSize );
478#endif
479        STICKYPROT( "(+" );
480        poly mon= nf[k];
481        while ( mon != NULL ) {
482            BOOLEAN found = FALSE;
483            int b= 0;
484            while ( found == FALSE ) {
485                if ( pEqual( mon, basis[b] ) )
486                    found= TRUE;
487                else
488                    b++;
489            }
490            number coeff = nCopy( pGetCoeff( mon ) );
491            v[k].setelem( b+1, coeff );
492            mon= pIter( mon );
493        }
494        STICKYPROT( ")" );
495    }
496    // gauss reduce
497    gaussReducer gauss( basisSize );
498    BOOLEAN isZero = FALSE;
499    fglmVector p;
500    for ( k= 0; (k < numMonoms) && (isZero == FALSE); k++ ) {
501        STICKYPROT( "(-" );
502        if ( ( isZero= gauss.reduce( v[k] )) == TRUE )
503            p= gauss.getDependence();
504        else
505            gauss.store();
506        STICKYPROT( ")" );
507    }
508    poly comb = NULL;
509    if ( isZero == TRUE ) {
510        number gcd = p.gcd();
511        if ( ! nIsZero( gcd ) && ! ( nIsOne( gcd ) ) )
512            p/= gcd;
513        nDelete( & gcd );
514        for ( k= 1; k <= p.size(); k++ ) {
515            if ( ! p.elemIsZero( k ) ) {
516                poly temp = pCopy( m[k-1] );
517                pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
518                comb= pAdd( comb, temp );
519            }
520        }
521        if ( ! nGreaterZero( pGetCoeff( comb ) ) ) comb= pNeg( comb );
522    }
523
524    // Free Memory
525    for ( k= 0; k < numMonoms; k++ ) {
526        pDelete( m + k );
527        pDelete( nf + k );
528    }
529    Free( (ADDRESS)m, numMonoms * sizeof( poly ) );
530    Free( (ADDRESS)nf, numMonoms * sizeof( poly ) );
531    // Warning: At this point all Vectors in v have to be initialized
532    for ( k= numMonoms - 1; k >= 0; k-- ) v[k].~fglmVector();
533    Free( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
534    for ( k= 0; k < basisSize; k++ )
535        pDelete( basis + k );
536    Free( (ADDRESS)basis, basisMax * sizeof( poly ) );
537    STICKYPROT( "\n" );
538    return comb;
539}
540
541#endif
542// Local Variables: ***
543// compile-command: "make Singular" ***
544// page-delimiter: "^\\(\\|//!\\)" ***
545// fold-internal-margins: nil ***
546// End: ***
Note: See TracBrowser for help on using the repository browser.