source: git/kernel/fglmzero.cc @ 0f401f

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