source: git/Singular/fglmcomb.cc @ d14712

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