source: git/kernel/fglmzero.cc @ fbc7cb

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