source: git/Singular/fglmzero.cc @ 17e692

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