source: git/Singular/fglmzero.cc @ ce9c2d

spielwiese
Last change on this file since ce9c2d was ce9c2d, checked in by Hans Schönemann <hannes@…>, 25 years ago
*hannes: aligne structures git-svn-id: file:///usr/local/Singular/svn/trunk@2874 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 34.3 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fglmzero.cc,v 1.19 1999-02-26 15:32:04 Singular Exp $
3
4/****************************************
5*  Computer Algebra System SINGULAR     *
6****************************************/
7/*
8* ABSTRACT - The FGLM-Algorithm
9*   Implementation of the fglm algorithm for 0-dimensional ideals,
10*   based on an idea by Faugere/Gianni/Lazard and Mora.
11*   The procedure CalculateFunctionals calculates the functionals
12*   which define the given ideal in the source ring. They build the
13*   input for GroebnerViaFunctionals, which defines the reduced
14*   groebner basis for the ideal in the destination ring.
15*/
16
17/* Changes:
18 * o FindUnivariatePolys added
19 */
20
21#include "mod2.h"
22
23#ifdef HAVE_FGLM
24#include "tok.h"
25#include "structs.h"
26#include "polys.h"
27#include "ideals.h"
28#include "ring.h"
29#include "ipid.h"
30#include "febase.h"
31#include "maps.h"
32#include "mmemory.h"
33#include "fglm.h"
34#include "fglmvec.h"
35#include "fglmgauss.h"
36// assumes, that NOSTREAMIO is set in factoryconf.h, which is included
37// by templates/list.h.
38#include <templates/ftmpl_list.h>
39#define PROT(msg)
40#define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
41#define PROT2(msg,arg)
42#define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
43#define fglmASSERT(ignore1,ignore2)
44
45
46// internal Version: 1.3.1.12
47// ============================================================
48//!      The idealFunctionals
49// ============================================================
50
51struct matElem
52{
53    int row;
54    number elem;
55};
56
57struct matHeader
58{
59    int size;
60    BOOLEAN owner;
61    matElem * elems;
62};
63
64class idealFunctionals
65{
66private:
67    int _block;
68    int _max;
69    int _size;
70    int _nfunc;
71    int * currentSize;
72    matHeader ** func;
73    matHeader * grow( int var );
74public:
75    idealFunctionals( int blockSize, int numFuncs );
76    ~idealFunctionals();
77
78    int dimen() const { fglmASSERT( _size>0, "called to early"); return _size; }
79    void endofConstruction();
80    void map( ring source );
81    void insertCols( int * divisors, int to );
82    void insertCols( int * divisors, const fglmVector to );
83    fglmVector addCols( const int var, int basisSize, const fglmVector v ) const;
84    fglmVector multiply( const fglmVector v, int var ) const;
85};
86
87idealFunctionals::idealFunctionals( int blockSize, int numFuncs )
88{
89    int k;
90    _block= blockSize;
91    _max= _block;
92    _size= 0;
93    _nfunc= numFuncs;
94
95    currentSize= (int *)Alloc( _nfunc*sizeof( int ) );
96    for ( k= _nfunc-1; k >= 0; k-- )
97        currentSize[k]= 0;
98
99    func= (matHeader **)Alloc( _nfunc*sizeof( matHeader * ) );
100    for ( k= _nfunc-1; k >= 0; k-- )
101        func[k]= (matHeader *)Alloc( _max*sizeof( matHeader ) );
102}
103
104idealFunctionals::~idealFunctionals()
105{
106    int k;
107    int l;
108    int row;
109    matHeader * colp;
110    matElem * elemp;
111    for ( k= _nfunc-1; k >= 0; k-- ) {
112        for ( l= _size-1, colp= func[k]; l >= 0; l--, colp++ ) {
113            if ( ( colp->owner == TRUE ) && ( colp->size > 0 ) ) {
114                for ( row= colp->size-1, elemp= colp->elems; row >= 0; row--, elemp++ )
115                    nDelete( & elemp->elem );
116                Free( (ADDRESS)colp->elems, colp->size*sizeof( matElem ) );
117            }
118        }
119        Free( (ADDRESS)func[k], _max*sizeof( matHeader ) );
120    }
121    Free( (ADDRESS)func, _nfunc*sizeof( matHeader * ) );
122    Free( (ADDRESS)currentSize, _nfunc*sizeof( int ) );
123}
124
125void
126idealFunctionals::endofConstruction()
127{
128    _size= currentSize[0];
129}
130
131void
132idealFunctionals::map( ring source )
133{
134    // maps from ring source to currentRing.
135    int var, col, row;
136    matHeader * colp;
137    matElem * elemp;
138    number newelem;
139
140    int * perm = (int *)Alloc0( (_nfunc+1)*sizeof( int ) );
141    maFindPerm( source->names, source->N, NULL, 0, currRing->names, currRing->N, NULL, 0, perm, NULL );
142    nSetMap( source->ch, source->parameter, source->P, source->minpoly );
143
144    matHeader ** temp = (matHeader **)Alloc( _nfunc*sizeof( matHeader * ));
145    for ( var= 0; var < _nfunc; var ++ ) {
146        for ( col= 0, colp= func[var]; col < _size; col++, colp++ ) {
147            if ( colp->owner == TRUE ) {
148                for ( row= colp->size-1, elemp= colp->elems; row >= 0; row--, elemp++ ) {
149                    newelem= nMap( elemp->elem );
150                    nDelete( & elemp->elem );
151                    elemp->elem= newelem;
152                }
153            }
154        }
155        temp[ perm[var+1]-1 ]= func[var];
156    }
157    Free( (ADDRESS)func, _nfunc*sizeof( matHeader * ) );
158    Free( (ADDRESS)perm, (_nfunc+1)*sizeof( int ) );
159    func= temp;
160}
161
162matHeader *
163idealFunctionals::grow( int var )
164{
165    if ( currentSize[var-1] == _max ) {
166        int k;
167        for ( k= _nfunc; k > 0; k-- )
168            func[k-1]= (matHeader *)ReAlloc( func[k-1], _max*sizeof( matHeader ), (_max + _block)*sizeof( matHeader ) );
169        _max+= _block;
170    }
171    currentSize[var-1]++;
172    return func[var-1] + currentSize[var-1] - 1;
173}
174
175void
176idealFunctionals::insertCols( int * divisors, int to )
177{
178    fglmASSERT( 0 < divisors[0] && divisors[0] <= _nfunc, "wrong number of divisors" );
179    int k;
180    BOOLEAN owner = TRUE;
181    matElem * elems = (matElem *)Alloc( sizeof( matElem ) );
182    elems->row= to;
183    elems->elem= nInit( 1 );
184    for ( k= divisors[0]; k > 0; k-- ) {
185        fglmASSERT( 0 < divisors[k] && divisors[k] <= _nfunc, "wrong divisor" );
186        matHeader * colp = grow( divisors[k] );
187        colp->size= 1;
188        colp->elems= elems;
189        colp->owner= owner;
190        owner= FALSE;
191    }
192}
193
194
195void
196idealFunctionals::insertCols( int * divisors, const fglmVector to )
197{
198    // divisors runs from divisors[0]..divisors[size-1]
199    fglmASSERT( 0 < divisors[0] && divisors[0] <= _nfunc, "wrong number of divisors" );
200    int k, l;
201    int numElems= to.numNonZeroElems();
202    matElem * elems;
203    matElem * elemp;
204    BOOLEAN owner = TRUE;
205    if ( numElems > 0 ) {
206        elems= (matElem *)Alloc( numElems * sizeof( matElem ) );
207        for ( k= 1, l= 1, elemp= elems; k <= numElems; k++, elemp++ ) {
208            while ( nIsZero( to.getconstelem(l) ) ) l++;
209            elemp->row= l;
210            elemp->elem= nCopy( to.getconstelem( l ) );
211            l++; // hochzaehlen, damit wir nicht noch einmal die gleiche Stelle testen
212        }
213    }
214    else
215        elems= NULL;
216    for ( k= divisors[0]; k > 0; k-- ) {
217        fglmASSERT( 0 < divisors[k] && divisors[k] <= _nfunc, "wrong divisor" );
218        matHeader * colp = grow( divisors[k] );
219        colp->size= numElems;
220        colp->elems= elems;
221        colp->owner= owner;
222        owner= FALSE;
223    }
224}
225
226fglmVector
227idealFunctionals::addCols( const int var, int basisSize, const fglmVector v ) const
228{
229    fglmVector result( basisSize );
230    matHeader * colp;
231    matElem * elemp;
232    number factor, temp;
233    int k, l;
234    int vsize = v.size();
235
236    fglmASSERT( currentSize[var-1]+1 >= vsize, "wrong v.size()" );
237    for ( k= 1, colp= func[var-1]; k <= vsize; k++, colp++ ) {
238        factor= v.getconstelem( k );
239        if ( ! nIsZero( factor ) ) {
240            for ( l= colp->size-1, elemp= colp->elems; l >= 0; l--, elemp++ ) {
241                temp= nMult( factor, elemp->elem );
242                number newelem= nAdd( result.getconstelem( elemp->row ), temp );
243                nDelete( & temp );
244                nNormalize( newelem );
245                result.setelem( elemp->row, newelem );
246            }
247        }
248    }
249    return result;
250}
251
252fglmVector
253idealFunctionals::multiply( const fglmVector v, int var ) const
254{
255    fglmASSERT( v.size() == _size, "multiply: v has wrong size");
256    fglmVector result( _size );
257    matHeader * colp;
258    matElem * elemp;
259    number factor, temp;
260    int k, l;
261    for ( k= 1, colp= func[var-1]; k <= _size; k++, colp++ ) {
262        factor= v.getconstelem( k );
263        if ( ! nIsZero( factor ) ) {
264            for ( l= colp->size-1, elemp= colp->elems; l >= 0; l--, elemp++ ) {
265                temp= nMult( factor, elemp->elem );
266                number newelem= nAdd( result.getconstelem( elemp->row ), temp );
267                nDelete( & temp );
268                nNormalize( newelem );
269                result.setelem( elemp->row, newelem );
270            }
271        }
272    }
273    return result;
274}
275
276// ============================================================
277//!      The old basis
278// ============================================================
279
280//     Contains the data for a border-monomial, i.e. the monom itself
281//      ( as a Singular-poly ) and its normal-form, described by an
282//      fglmVector. The size of the Vector depends on the current size of
283//      the basis when the normalForm was computed.
284//     monom gets deleted when borderElem comes out of scope.
285class borderElem
286{
287public:
288    poly monom;
289    fglmVector nf;
290    borderElem() : monom(NULL), nf() {}
291    borderElem( poly p, fglmVector n ) : monom( p ), nf( n ) {}
292    ~borderElem() { pDelete1( &monom ); }
293#ifndef HAVE_EXPLICIT_CONSTR
294    void insertElem( poly p, fglmVector n )
295    {
296        monom= p;
297        nf= n;
298    }
299#endif
300};
301
302//     This class contains the relevant data for the 'candidates'
303//     The declaration of class fglmSelem is found in fglm.h
304
305fglmSelem::fglmSelem( poly p, int var ) : monom( p ), numVars( 0 )
306{
307    for ( int k = pVariables; k > 0; k-- )
308        if ( pGetExp( monom, k ) > 0 )
309            numVars++;
310    divisors= (int *)Alloc( (numVars+1)*sizeof( int ) );
311    divisors[0]= 0;
312    newDivisor( var );
313}
314
315void
316fglmSelem::cleanup()
317{
318    Free( (ADDRESS)divisors, (numVars+1)*sizeof( int ) );
319}
320
321//     The data-structure for the Functional-Finding-Algorithm.
322class fglmSdata
323{
324private:
325    ideal theIdeal;
326    int idelems;
327
328    int basisBS;
329    int basisMax;
330    int basisSize;
331    polyset basis;  //. rem: runs from basis[1]..basis[dimen]
332
333    int borderBS;
334    int borderMax;
335    int borderSize;
336    borderElem * border;  //. rem: runs from border[1]..border[dimen]
337
338    List<fglmSelem> nlist;
339    BOOLEAN _state;
340public:
341    fglmSdata( const ideal thisIdeal );
342    ~fglmSdata();
343
344    BOOLEAN state() const { return _state; };
345    int getBasisSize() const { return basisSize; };
346    int newBasisElem( poly & p );
347    void newBorderElem( poly & m, fglmVector v );
348    BOOLEAN candidatesLeft() const { return ( nlist.isEmpty() ? FALSE : TRUE ); }
349    fglmSelem nextCandidate();
350    void updateCandidates();
351    int getEdgeNumber( const poly m ) const;
352    poly getSpanPoly( int number ) const { return pCopy( (theIdeal->m)[number-1] ); }
353    fglmVector getVectorRep( const poly m );
354    fglmVector getBorderDiv( const poly m, int & var ) const;
355};
356
357fglmSdata::fglmSdata( const ideal thisIdeal )
358{
359    // An dieser Stelle kann die BlockSize ( =BS ) noch sinnvoller berechnet
360    // werden, jenachdem wie das Ideal aussieht.
361    theIdeal= thisIdeal;
362    idelems= IDELEMS( theIdeal );
363
364    basisBS= 100;
365    basisMax= basisBS;
366    basisSize= 0;
367    basis= (polyset)Alloc( basisMax*sizeof( poly ) );
368
369    borderBS= 100;
370    borderMax= borderBS;
371    borderSize= 0;
372#ifndef HAVE_EXPLICIT_CONSTR
373    border= new borderElem[ borderMax ];
374#else
375    border= (borderElem *)Alloc( borderMax*sizeof( borderElem ) );
376#endif
377    // rem: the constructors are called in newBorderElem().
378    _state= TRUE;
379}
380
381fglmSdata::~fglmSdata()
382{
383    for ( int k = basisSize; k > 0; k-- )
384        pDelete1( basis + k );  //. rem: basis runs from basis[1]..basis[basisSize]
385    Free( (ADDRESS)basis, basisMax*sizeof( poly ) );
386#ifndef HAVE_EXPLICIT_CONSTR
387    delete [] border;
388#else
389    for ( int l = borderSize; l > 0; l-- )
390        // rem: the polys of borderElem are deleted via ~borderElem()
391        border[l].~borderElem();
392    Free( (ADDRESS)border, borderMax*sizeof( borderElem ) );
393#endif
394}
395
396//     Inserts poly p without copying into basis, reAllocs Memory if necessary,
397//      increases basisSize and returns the new basisSize.
398//     Remember: The elements of basis are deleted via pDelete in ~fglmSdata!
399//     Sets m= NULL to indicate that now basis is ow(e?)ing the poly.
400int
401fglmSdata::newBasisElem( poly & m )
402{
403    basisSize++;
404    if ( basisSize == basisMax ) {
405        basis= (polyset)ReAlloc( basis, basisMax*sizeof( poly ), (basisMax + basisBS)*sizeof( poly ) );
406        basisMax+= basisBS;
407    }
408    basis[basisSize]= m;
409    m= NULL;
410    return basisSize;
411}
412
413//     Inserts poly p and fglmvector v without copying into border, reAllocs Memory
414//      if necessary, and increases borderSize.
415//     Remember: The poly of border is deleted via ~borderElem in ~fglmSdata!
416//     Sets m= NULL to indicate that now border is ow(e?)ing the poly.
417void
418fglmSdata::newBorderElem( poly & m, fglmVector v )
419{
420    borderSize++;
421    if ( borderSize == borderMax ) {
422#ifndef HAVE_EXPLICIT_CONSTR
423        borderElem * tempborder = new borderElem[ borderMax+borderBS ];
424        for ( int k = 0; k < borderMax; k++ ) {
425            tempborder[k]= border[k];
426            border[k].insertElem( NULL, fglmVector() );
427        }
428        delete [] border;
429        border= tempborder;
430#else
431        border= (borderElem *)ReAlloc( border, borderMax*sizeof( borderElem ), (borderMax + borderBS)*sizeof( borderElem ) );
432#endif
433        borderMax+= borderBS;
434    }
435#ifndef HAVE_EXPLICIT_CONSTR
436    border[borderSize].insertElem( m, v );
437#else
438    border[borderSize].borderElem( m, v );
439#endif
440    m= NULL;
441}
442
443fglmSelem
444fglmSdata::nextCandidate()
445{
446    fglmSelem result = nlist.getFirst();
447    nlist.removeFirst();
448    return result;
449}
450
451//     Multiplies basis[basisSize] with all ringvariables and inserts the new monomials
452//      into the list of candidates, according to the given order. If a monomial already
453//      exists, then "insertions" and "divisors" are updated.
454//     Assumes that ringvar(k) < ringvar(l) for k > l
455void
456fglmSdata::updateCandidates()
457{
458    ListIterator<fglmSelem> list = nlist;
459    fglmASSERT( basisSize > 0 && basisSize < basisMax, "Error(1) in fglmSdata::updateCandidates - wrong bassSize" );
460    poly m = basis[basisSize];
461    poly newmonom = NULL;
462    int k = pVariables;
463    BOOLEAN done = FALSE;
464    int state = 0;
465    while ( k >= 1 ) {
466        newmonom = pCopy( m );
467        pIncrExp( newmonom, k );
468        pSetm( newmonom );
469        done= FALSE;
470        while ( list.hasItem() && (done == FALSE) ) {
471            if ( (state= pComp( list.getItem().monom, newmonom )) < 0 )
472                list++;
473            else done= TRUE;
474        }
475        if ( done == FALSE ) {
476            nlist.append( fglmSelem( newmonom, k ) );
477            break;
478        }
479        if ( state == 0 ) {
480            list.getItem().newDivisor( k );
481            pDelete1( &newmonom );
482        }
483        else {
484            list.insert( fglmSelem( newmonom, k ) );
485        }
486        k--;
487    }
488    while ( --k >= 1 ) {
489        newmonom= pCopy( m ); // HIER
490        pIncrExp( newmonom, k );
491        pSetm( newmonom );
492        nlist.append( fglmSelem( newmonom, k ) );
493    }
494}
495
496//     if p == pHead( (theIdeal->m)[k] ) return k, 0 otherwise
497//     (Assumes that pEqual just checks the leading monomials without
498//      coefficients.)
499int
500fglmSdata::getEdgeNumber( const poly m ) const
501{
502    for ( int k = idelems; k > 0; k-- )
503        if ( pEqual( m, (theIdeal->m)[k-1] ) )
504            return k;
505    return 0;
506}
507
508//     Returns the fglmVector v, s.t.
509//        p = v[1]*basis(1) + .. + v[basisSize]*basis(basisSize)
510//     So the size of v depends on the current size of the basis.
511//     Assumes that such an representation exists, i.e. all monoms in p have to be
512//      smaller than basis[basisSize] and that basis[k] < basis[l] for k < l.
513fglmVector
514fglmSdata::getVectorRep( const poly p )
515{
516    fglmVector temp( basisSize );
517    poly m = p;
518    int num = basisSize;
519    while ( m != NULL ) {
520        int comp = pComp( m, basis[num] );
521        if ( comp == 0 ) {
522            fglmASSERT( num > 0, "Error(1) in fglmSdata::getVectorRep" );
523            number newelem = nCopy( pGetCoeff( m ) );
524            temp.setelem( num, newelem );
525            num--;
526            pIter( m );
527        }
528        else {
529            if ( comp < 0 ) {
530                num--;
531            }
532            else {
533                // This is the place where we can detect if the sourceIdeal
534                // is not reduced. In this case m is not in basis[]. Since basis[]
535                // is ordered this is then and only then the case, if basis[i]<m
536                // and basis[j]>m for all j>i
537                _state= FALSE;
538                return temp;
539            }
540        }
541    }
542    return temp;
543}
544
545//     Searches through the border for a monomoial bm which devides m and returns
546//      its normalform in vector representation.
547//     var contains the number of the variable v, s.t. bm = m * v
548fglmVector
549fglmSdata::getBorderDiv( const poly m, int & var ) const
550{
551//     int num2 = borderSize;
552//     while ( num2 > 0 ) {
553//         poly temp = border[num2].monom;
554//         if ( pDivisibleBy( temp, m ) ) {
555//             poly divisor = pDivideM( m, temp );
556//             int var = pIsPurePower( divisor );
557//             if ( (var != 0) && (pGetCoeff( divisor, var ) == 1) ) {
558//                 Print( "poly %s divides poly %s", pString( temp ), pString( m ) );
559//             }
560//         }
561//         num2--;
562//     }
563    int num = borderSize;
564    while ( num > 0 ) {
565        poly temp = border[num].monom;
566        if ( pDivisibleBy( temp, m ) ) {
567            var = pVariables;
568            while ( var > 0 ) {
569                if ( (pGetExp( m, var ) - pGetExp( temp, var )) == 1 )
570                    return border[num].nf;
571                var--;
572            }
573        }
574        num--;
575    }
576    return fglmVector();
577}
578
579//     Calculates the defining Functionals for the ideal "theIdeal" and
580//     returns them in "l".
581//     The ideal has to be zero-dimensional and reduced and has to be a
582//     real subset of the polynomal ring.
583//     In any case it has to be zero-dimensional and minimal (check this
584//      via fglmIdealcheck). Any minimal but not reduced ideal is detected.
585//      In this case it returns FglmNotReduced.
586//     If the base domain is Q, the leading coefficients of the polys
587//     have to be in Z.
588//     returns TRUE if the result is valid, FALSE if theIdeal
589//      is not reduced.
590static BOOLEAN
591CalculateFunctionals( const ideal & theIdeal, idealFunctionals & l )
592{
593    fglmSdata data( theIdeal );
594
595    // insert pOne() into basis and update the workingList:
596    poly one = pOne();
597    data.newBasisElem( one );
598    data.updateCandidates();
599
600    STICKYPROT(".");
601    while ( data.candidatesLeft() == TRUE ) {
602        fglmSelem candidate = data.nextCandidate();
603        if ( candidate.isBasisOrEdge() == TRUE ) {
604            int edge = data.getEdgeNumber( candidate.monom );
605            if ( edge != 0 ) {
606                // now candidate is an edge, i.e. we know its normalform:
607                // NF(p) = - ( tail(p)/LC(p) )
608                poly nf = data.getSpanPoly( edge );
609                pNorm( nf );
610                pDelete1( &nf );  //. deletes the leadingmonomial
611                nf= pNeg( nf );
612                fglmVector nfv = data.getVectorRep( nf );
613                l.insertCols( candidate.divisors, nfv );
614                data.newBorderElem( candidate.monom, nfv );
615                pDelete( &nf );
616                STICKYPROT( "+" );
617            }
618            else {
619                int basis= data.newBasisElem( candidate.monom );
620                data.updateCandidates();
621                l.insertCols( candidate.divisors, basis );
622                STICKYPROT( "." );
623            }
624        }
625        else {
626            int var = 0;
627            fglmVector temp = data.getBorderDiv( candidate.monom, var );
628            fglmASSERT( var > 0, "this should never happen" );
629            fglmVector nfv = l.addCols( var, data.getBasisSize(), temp );
630            data.newBorderElem( candidate.monom, nfv );
631            l.insertCols( candidate.divisors, nfv );
632            STICKYPROT( "-" );
633        }
634        candidate.cleanup();
635    } //. while ( data.candidatesLeft() == TRUE )
636    l.endofConstruction();
637    STICKYPROT2( "\nvdim= %i\n", data.getBasisSize() );
638    return ( data.state() );
639}
640
641// ============================================================
642//!      The new basis
643// ============================================================
644
645//     The declaration of class fglmDelem is found in fglm.h
646
647fglmDelem::fglmDelem( poly & m, fglmVector mv, int v ) : v( mv ), insertions( 0 ), var( v )
648{
649    monom= m;
650    m= NULL;
651    for ( int k = pVariables; k > 0; k-- )
652        if ( pGetExp( monom, k ) > 0 )
653            insertions++;
654    // Wir gehen davon aus, dass ein fglmDelem direkt bei der Erzeugung
655    // auch in eine Liste eingefuegt wird. Daher wird hier automatisch
656    // newDivisor aufgerufen ( v teilt ja m )
657    newDivisor();
658}
659
660void
661fglmDelem::cleanup()
662{
663    if ( monom != NULL ) {
664        pDelete1( &monom );
665    }
666}
667
668class oldGaussElem
669{
670public:
671    fglmVector v;
672    fglmVector p;
673    number pdenom;
674    number fac;
675
676#ifndef HAVE_EXPLICIT_CONSTR
677    oldGaussElem() : v(), p(), pdenom( NULL ), fac( NULL ) {}
678#endif
679    oldGaussElem( const fglmVector newv, const fglmVector newp, number & newpdenom, number & newfac ) : v( newv ), p( newp ), pdenom( newpdenom ), fac( newfac )
680    {
681        newpdenom= NULL;
682        newfac= NULL;
683    }
684    ~oldGaussElem();
685#ifndef HAVE_EXPLICIT_CONSTR
686    void insertElem( const fglmVector newv, const fglmVector newp, number & newpdenom, number & newfac )
687    {
688        v= newv;
689        p= newp;
690        pdenom= newpdenom;
691        fac= newfac;
692        newpdenom= NULL;
693        newfac= NULL;
694    }
695#endif
696};
697
698oldGaussElem::~oldGaussElem()
699{
700    nDelete( & fac );
701    nDelete( & pdenom );
702}
703
704
705class fglmDdata
706{
707private:
708    int dimen;
709    oldGaussElem * gauss;
710    BOOLEAN * isPivot;  // [1]..[dimen]
711    int * perm;  // [1]..[dimen]
712    int basisSize;  //. the CURRENT basisSize, i.e. basisSize <= dimen
713    polyset basis;  // [1]..[dimen]. The monoms of the new Vectorspace-basis
714
715    int groebnerBS;
716    int groebnerSize;
717    ideal destId;
718
719    List<fglmDelem> nlist;
720public:
721    fglmDdata( int dimension );
722    ~fglmDdata();
723
724    int getBasisSize() const { return basisSize; }
725    BOOLEAN candidatesLeft() const { return ( nlist.isEmpty() ? FALSE : TRUE ); }
726    fglmDelem nextCandidate();
727    void newBasisElem( poly & m, fglmVector v, fglmVector p, number & denom );
728    void updateCandidates( poly m, const fglmVector v );
729    void newGroebnerPoly( fglmVector & v, poly & p );
730    void gaussreduce( fglmVector & v, fglmVector & p, number & denom );
731    ideal buildIdeal()
732    {
733        idSkipZeroes( destId );
734        return destId;
735    }
736};
737
738fglmDdata::fglmDdata( int dimension )
739{
740    int k;
741    dimen= dimension;
742    basisSize= 0;
743    //. All arrays run from [1]..[dimen], thus Alloc( dimen + 1 )!
744#ifndef HAVE_EXPLICIT_CONSTR
745    gauss= new oldGaussElem[ dimen+1 ];
746#else
747    gauss= (oldGaussElem *)Alloc( (dimen+1)*sizeof( oldGaussElem ) );
748#endif
749    isPivot= (BOOLEAN *)Alloc( (dimen+1)*sizeof( BOOLEAN ) );
750    for ( k= dimen; k > 0; k-- ) isPivot[k]= FALSE;
751    perm= (int *)Alloc( (dimen+1)*sizeof( int ) );
752    basis= (polyset)Alloc( (dimen+1)*sizeof( poly ) );
753    groebnerBS= 16;
754    groebnerSize= 0;
755    destId= idInit( groebnerBS, 1 );
756}
757
758fglmDdata::~fglmDdata()
759{
760    fglmASSERT( dimen == basisSize, "Es wurden nicht alle BasisElemente gefunden!" );
761    int k;
762#ifndef HAVE_EXPLICIT_CONSTR
763    delete [] gauss;
764#else
765    for ( k= dimen; k > 0; k-- )
766        gauss[k].~oldGaussElem();
767    Free( (ADDRESS)gauss, (dimen+1)*sizeof( oldGaussElem ) );
768#endif
769    Free( (ADDRESS)isPivot, (dimen+1)*sizeof( BOOLEAN ) );
770    Free( (ADDRESS)perm, (dimen+1)*sizeof( int ) );
771    //. Remember: There is no poly in basis[0], thus k > 0
772    for ( k= dimen; k > 0; k-- )
773        pDelete1( basis + k );
774    Free( (ADDRESS)basis, (dimen+1)*sizeof( poly ) );
775}
776
777fglmDelem
778fglmDdata::nextCandidate()
779{
780    fglmDelem result = nlist.getFirst();
781    nlist.removeFirst();
782    return result;
783}
784
785void
786fglmDdata::newBasisElem( poly & m, fglmVector v, fglmVector p, number & denom )
787{
788// inserts m as a new basis monom. m is NOT copied but directly inserted.
789// returns m=NULL to indicate, that now basis is oweing m.
790    basisSize++;
791    basis[basisSize]= m;
792    m= NULL;
793    int k= 1;
794    while ( nIsZero(v.getconstelem(k)) || isPivot[k] ) {
795        k++;
796    }
797    fglmASSERT( k <= dimen, "Error(1) in fglmDdata::pivot-search");
798    number pivot= v.getconstelem( k );
799    int pivotcol = k;
800    k++;
801    while ( k <= dimen ) {
802        if ( ! nIsZero( v.getconstelem(k) ) && ! isPivot[k] ) {
803            if ( nGreater( v.getconstelem( k ), pivot ) ) {
804                pivot= v.getconstelem( k );
805                pivotcol= k;
806            }
807        }
808        k++;
809    }
810    fglmASSERT( ! nIsZero( pivot ), "Error(2) fglmDdata::Pivotelement ist Null" );
811    isPivot[ pivotcol ]= TRUE;
812    perm[basisSize]= pivotcol;
813
814    pivot= nCopy( v.getconstelem( pivotcol ) );
815#ifndef HAVE_EXPLICIT_CONSTR
816    gauss[basisSize].insertElem( v, p, denom, pivot );
817#else
818    gauss[basisSize].oldGaussElem( v, p, denom, pivot );
819#endif
820}
821
822void
823fglmDdata::updateCandidates( poly m, const fglmVector v )
824{
825    ListIterator<fglmDelem> list = nlist;
826    poly newmonom = NULL;
827    int k = pVariables;
828    BOOLEAN done = FALSE;
829    int state = 0;
830    while ( k >= 1 ) {
831        newmonom = pCopy( m );
832        pIncrExp( newmonom, k );
833        pSetm( newmonom );
834        done= FALSE;
835        while ( list.hasItem() && (done == FALSE) ) {
836            if ( (state= pComp( list.getItem().monom, newmonom )) < 0 )
837                list++;
838            else done= TRUE;
839        }
840        if ( done == FALSE ) {
841            nlist.append( fglmDelem( newmonom, v, k ) );
842            break;
843        }
844        if ( state == 0 ) {
845            list.getItem().newDivisor();
846            pDelete1( & newmonom );
847        }
848        else {
849            list.insert( fglmDelem( newmonom, v, k ) );
850        }
851        k--;
852    }
853    while ( --k >= 1 ) {
854        newmonom= pCopy( m );
855        pIncrExp( newmonom, k );
856        pSetm( newmonom );
857        nlist.append( fglmDelem( newmonom, v, k ) );
858    }
859}
860
861void
862fglmDdata::newGroebnerPoly( fglmVector & p, poly & m )
863// Inserts gp = p[1]*basis(1)+..+p[basisSize]*basis(basisSize)+p[basisSize+1]*m as
864//  a new groebner polynomial for the ideal.
865// All elements (monomials and coefficients) of gp are copied, instead of m.
866// Assumes that p.length() == basisSize+1.
867{
868    //. Baue das Polynom von oben nach unten:
869    fglmASSERT( p.size() == basisSize+1, "GP::newGroebnerPoly: p has wrong size" );
870    int k;
871    poly result = m;
872    poly temp = result;
873    m= NULL;
874    if ( nGetChar() > 0 ) {
875        number lead = nCopy( p.getconstelem( basisSize+1 ) );
876        p /= lead;
877        nDelete( & lead );
878    }
879    if ( nGetChar() == 0 ) {
880        number gcd= p.gcd();
881        fglmASSERT( ! nIsZero( gcd ), "FATAL: gcd and thus p is zero" );
882        if ( ! nIsOne( gcd ) )
883            p /= gcd;
884        nDelete( & gcd );
885    }
886    pSetCoeff( result, nCopy( p.getconstelem( basisSize+1 ) ) );
887    for ( k= basisSize; k > 0; k-- ) {
888        if ( ! nIsZero( p.getconstelem( k ) ) ) {
889            temp->next= pCopy( basis[k] );
890            pIter( temp );
891            pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
892        }
893    }
894    pSetm( result );
895    if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
896    if ( groebnerSize == IDELEMS( destId ) ) {
897        pEnlargeSet( & destId->m, IDELEMS( destId ), groebnerBS );
898        IDELEMS( destId )+= groebnerBS;
899    }
900    (destId->m)[groebnerSize]= result;
901    groebnerSize++;
902}
903
904void
905fglmDdata::gaussreduce( fglmVector & v, fglmVector & p, number & pdenom )
906{
907    int k;
908    number fac1, fac2;
909    number temp;
910    fglmASSERT( pdenom == NULL, "pdenom in gaussreduce should be NULL" );
911    pdenom= nInit( 1 );
912    number vdenom = v.clearDenom();
913    if ( ! nIsOne( vdenom ) && ! nIsZero( vdenom ) ) {
914        p.setelem( p.size(), vdenom );
915    }
916    else {
917        nDelete( &vdenom );
918    }
919    number gcd = v.gcd();
920    if ( ! nIsOne( gcd ) && ! nIsZero( gcd ) ) {
921        v /= gcd;
922        number temp= nMult( pdenom, gcd );
923        nDelete( &pdenom );
924        pdenom= temp;
925    }
926    nDelete( & gcd );
927
928    for ( k= 1; k <= basisSize; k++ ) {
929
930        if ( ! v.elemIsZero( perm[k] ) ) {
931            fac1= gauss[k].fac;
932            fac2= nCopy( v.getconstelem( perm[k] ) );
933            v.nihilate( fac1, fac2, gauss[k].v );
934            fac1= nMult( fac1, gauss[k].pdenom );
935            temp= nMult( fac2, pdenom );
936            nDelete( &fac2 );
937            fac2= temp;
938              p.nihilate( fac1, fac2, gauss[k].p );
939            temp= nMult( pdenom, gauss[k].pdenom );
940            nDelete( &pdenom );
941            pdenom= temp;
942
943            nDelete( & fac1 );
944            nDelete( & fac2 );
945            number gcd = v.gcd();
946            if ( ! nIsOne( gcd ) && ! nIsZero( gcd ) ) {
947                v /= gcd;
948                number temp= nMult( pdenom, gcd );
949                nDelete( &pdenom );
950                pdenom= temp;
951            }
952            nDelete( & gcd );
953            gcd= p.gcd();
954            temp= nGcd( pdenom, gcd );
955            nDelete( &gcd );
956            gcd= temp;
957            if ( ! nIsZero( gcd ) && ! nIsOne( gcd ) ) {
958                p /= gcd;
959                temp= nDiv( pdenom, gcd );
960                nDelete( & pdenom );
961                pdenom= temp;
962                nNormalize( pdenom );
963            }
964            nDelete( & gcd );
965        }
966    }
967}
968
969static ideal
970GroebnerViaFunctionals( const idealFunctionals & l )
971// Calculates the groebnerBasis for the ideal which is defined by l.
972// The dimension of l has to be finite.
973// The result is in reduced form.
974{
975    fglmDdata data( l.dimen() );
976
977    // insert pOne() and update workinglist:
978    poly one = pOne();
979    data.updateCandidates( one, fglmVector(l.dimen(), 1) );
980    number nOne = nInit( 1 );
981    data.newBasisElem( one, fglmVector( l.dimen(), 1 ), fglmVector( 1, 1 ), nOne );
982    STICKYPROT( "." );
983    while ( data.candidatesLeft() == TRUE ) {
984        fglmDelem candidate = data.nextCandidate();
985        if ( candidate.isBasisOrEdge() == TRUE ) {
986            // Now we have the chance to find a new groebner polynomial
987
988            // v is the vector-representation of candidate.monom
989            // some elements of v are zeroed in data.gaussreduce(). Which
990            // ones and how this was done is stored in p.
991            // originalV containes the unchanged v, which is later inserted
992            // into the working list (via data.updateCandidates().
993            fglmVector v = l.multiply( candidate.v, candidate.var );
994            fglmVector originalV = v;
995            fglmVector p( data.getBasisSize()+1, data.getBasisSize()+1 );
996            number pdenom = NULL;
997            data.gaussreduce( v, p, pdenom );
998            if ( v.isZero() ) {
999                // Now v is linear dependend to the already found basis elements.
1000                // This means that v (rsp. candidate.monom) is the leading
1001                // monomial of the next groebner-basis polynomial.
1002                data.newGroebnerPoly( p, candidate.monom );
1003                nDelete( & pdenom );
1004                STICKYPROT( "+" );
1005            }
1006            else {
1007                // no linear dependence could be found, so v ( rsp. monom )
1008                // is a basis monomial. We store the zeroed version ( i.e. v
1009                // and not originalV ) as well as p, the denomiator and all
1010                // the other stuff.
1011                // erst updateCandidates, dann newBasisELem!!!
1012                data.updateCandidates( candidate.monom, originalV );
1013                data.newBasisElem( candidate.monom, v, p, pdenom );
1014                STICKYPROT( "." );
1015            }
1016        }
1017        else {
1018            STICKYPROT( "-" );
1019            candidate.cleanup();
1020        }
1021    }  //. while data.candidatesLeft()
1022    STICKYPROT( "\n" );
1023    return ( data.buildIdeal() );
1024}
1025//<-
1026
1027static ideal
1028FindUnivariatePolys( const idealFunctionals & l )
1029{
1030    fglmVector v;
1031    fglmVector p;
1032    ideal destIdeal = idInit( pVariables, 1 );
1033
1034    int i;
1035    BOOLEAN isZero;
1036    for ( i= 1; i <= pVariables; i++ ) {
1037        // main loop
1038        STICKYPROT2( "(%i)", i );
1039        gaussReducer gauss( l.dimen() );
1040        isZero= FALSE;
1041        v= fglmVector( l.dimen(), 1 );
1042        while ( isZero == FALSE ) {
1043            if ( (isZero= gauss.reduce( v )) == TRUE ) {
1044                STICKYPROT( "+" );
1045                p= gauss.getDependence();
1046                number gcd= p.gcd();
1047                if ( ! nIsOne( gcd ) ) {
1048                    p /= gcd;
1049                }
1050                nDelete( & gcd );
1051                int k;
1052                poly temp = NULL;
1053                poly result;
1054                for ( k= p.size(); k > 0; k-- ) {
1055                    number n = nCopy( p.getconstelem( k ) );
1056                    if ( ! nIsZero( n ) ) {
1057                        if ( temp == NULL ) {
1058                            result= pOne();
1059                            temp= result;
1060                        }
1061                        else {
1062                            temp->next= pOne();
1063                            pIter( temp );
1064                        }
1065                        pSetCoeff( temp, n );
1066                        pSetExp( temp, i, k-1 );
1067                        pSetm( temp );
1068                    }
1069                }
1070                if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
1071                (destIdeal->m)[i-1]= result;
1072            }
1073            else {
1074                STICKYPROT( "." );
1075                gauss.store();
1076                v= l.multiply( v, i );
1077            }
1078        }
1079    }
1080    STICKYPROT( "\n" );
1081    return destIdeal;
1082}
1083
1084// for a descritption of the parameters see fglm.h
1085BOOLEAN
1086fglmzero( idhdl sourceRingHdl, ideal & sourceIdeal, idhdl destRingHdl, ideal & destIdeal, BOOLEAN switchBack, BOOLEAN deleteIdeal )
1087{
1088    idhdl initialRingHdl = currRingHdl;
1089    BOOLEAN fglmok;
1090
1091    if ( currRingHdl != sourceRingHdl )
1092        rSetHdl( sourceRingHdl, TRUE );
1093    idealFunctionals L( 100, pVariables );
1094    fglmok = CalculateFunctionals( sourceIdeal, L );
1095    if ( deleteIdeal == TRUE )
1096        idDelete( & sourceIdeal );
1097    rSetHdl( destRingHdl, TRUE );
1098    if ( fglmok == TRUE ) {
1099        L.map( IDRING( sourceRingHdl ) );
1100        destIdeal= GroebnerViaFunctionals( L );
1101    }
1102    if ( (switchBack == TRUE) && (currRingHdl != initialRingHdl) )
1103        rSetHdl( initialRingHdl, TRUE );
1104    return fglmok;
1105}
1106
1107BOOLEAN
1108FindUnivariateWrapper( ideal source, ideal & destIdeal )
1109{
1110    BOOLEAN fglmok;
1111
1112    idealFunctionals L( 100, pVariables );
1113    fglmok = CalculateFunctionals( source, L );
1114    if ( fglmok == TRUE ) {
1115        destIdeal= FindUnivariatePolys( L );
1116        return TRUE;
1117    }
1118    else
1119        return FALSE;
1120}
1121
1122
1123#endif
1124
1125// ----------------------------------------------------------------------------
1126// Local Variables: ***
1127// compile-command: "make Singular" ***
1128// page-delimiter: "^\\(\\|//!\\)" ***
1129// fold-internal-margins: nil ***
1130// End: ***
Note: See TracBrowser for help on using the repository browser.