source: git/Singular/fglmzero.cc @ 4508ce5

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