source: git/Singular/fglmcomb.cc @ 7df4ee

spielwiese
Last change on this file since 7df4ee was 7df4ee, checked in by Hans Schönemann <hannes@…>, 26 years ago
*hannes: POSIX-fixes part2 git-svn-id: file:///usr/local/Singular/svn/trunk@2607 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.14 1998-10-28 12:43: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#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        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                pIter( *pptr );
63                temp= *pptr;
64            }
65            else {
66                todelete= temp;
67                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            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                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        pIter( temp );
198    }
199
200    nf= (polyset)Alloc( numMonoms * sizeof( poly ) );
201
202#ifndef HAVE_EXPLICIT_CONSTR
203    mv= new fglmVector[ numMonoms ];
204#else
205    mv= (fglmVector *)Alloc( numMonoms * sizeof( fglmVector ) );
206#endif
207
208#ifndef HAVE_EXPLICIT_CONSTR
209    v= new fglmVector[ numMonoms ];
210#else
211    v= (fglmVector *)Alloc( numMonoms * sizeof( fglmVector ) );
212#endif
213
214    basisMax= basisBS;
215    basis= (polyset)Alloc( basisMax * sizeof( poly ) );
216
217    weights= (int *)Alloc( IDELEMS( source ) * sizeof( int ) );
218    STICKYPROT( "weights: " );
219    for ( k= 0; k < IDELEMS( source ); k++ ) {
220        poly temp= (source->m)[k];
221        int w = 0;
222        while ( temp != NULL ) {
223            w+= nSize( pGetCoeff( temp ) );
224            pIter( temp );
225        }
226        weights[k]= w;
227        STICKYPROT2( "%i ", w );
228    }
229    STICKYPROT( "\n" );
230    lengthes= (int *)Alloc( numMonoms * sizeof( int ) );
231    order= (int *)Alloc( numMonoms * sizeof( int ) );
232
233    // calculate the NormalForm in a special way
234    for ( k= 0; k < numMonoms; k++ )
235    {
236        STICKYPROT( "#" );
237        poly current = pCopy( m[k] );
238        fglmVector currV( numMonoms, k+1 );
239
240        fglmReduce( & current, currV, m, numMonoms, source, weights );
241        poly temp = current;
242        int b;
243        while ( temp != NULL )
244        {
245            BOOLEAN found = FALSE;
246            for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
247            {
248                if ( pEqual( temp, basis[b] ) )
249                {
250                    found= TRUE;
251                }
252            }
253            if ( found == FALSE )
254            {
255                if ( basisSize == basisMax )
256                {
257                    // Expand the basis
258                    basis= (polyset)ReAlloc( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
259                    basisMax+= basisBS;
260                }
261//                 basis[basisSize]= pOne();
262//                 pSetExpV( basis[basisSize], temp->exp );
263//                 pSetm( basis[basisSize] );
264                basis[basisSize]= pNew();
265                pCopy2( basis[basisSize], 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 ( pEqual( 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    Free( (ADDRESS)nf, numMonoms * sizeof( poly ) );
308    Free( (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    Free( (ADDRESS)lengthes, numMonoms * sizeof( int ) );
399    Free( (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    Free( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
406#endif
407
408    for ( k= 0; k < basisSize; k++ )
409        pDelete( basis + k );
410    Free( (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    Free( (ADDRESS)mv, numMonoms * sizeof( fglmVector ) );
418#endif
419
420    for ( k= 0; k < numMonoms; k++ )
421        pDelete( m + k );
422    Free( (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)Alloc( numMonoms * sizeof( poly ) );
452    nf= (polyset)Alloc( numMonoms * sizeof( poly ) );
453    // Warning: The fglmVectors in v are yet not initialized
454    v= (fglmVector *)Alloc( numMonoms * sizeof( fglmVector ) );
455    basisMax= basisBS;
456    basis= (polyset)Alloc( 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 ( pEqual( sm, basis[b] ) ) found= TRUE;
478            if ( found == FALSE ) {
479                // Expand the basis
480                if ( basisSize == basisMax ) {
481                    basis= (polyset)ReAlloc( 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 ( pEqual( 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    Free( (ADDRESS)m, numMonoms * sizeof( poly ) );
554    Free( (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    Free( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
558    for ( k= 0; k < basisSize; k++ )
559        pDelete( basis + k );
560    Free( (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.