source: git/Singular/fglmzero.cc @ 5480da

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