source: git/Singular/fglmzero.cc @ 138f0c

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