source: git/kernel/fglmzero.cc @ fcb8022

spielwiese
Last change on this file since fcb8022 was e710e2, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: code cleanup git-svn-id: file:///usr/local/Singular/svn/trunk@7959 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.5 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: fglmzero.cc,v 1.3 2005-04-29 17:35:17 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 "structs.h"
25#include "polys.h"
26#include "ideals.h"
27#include "ring.h"
28#include "febase.h"
29#include "maps.h"
30#include "omalloc.h"
31#include "kstd1.h" // for kNF (see fglmquot)
32#include "intvec.h"
33#include "fglm.h"
34#include "fglmvec.h"
35#include "fglmgauss.h"
36// assumes, that NOSTREAMIO is set in factoryconf.h, which is included
37// by templates/list.h.
38#include <templates/ftmpl_list.h>
39#define PROT(msg)
40#define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg)
41#define PROT2(msg,arg)
42#define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg)
43#define fglmASSERT(ignore1,ignore2)
44
45
46
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() { pDeleteLm(&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        pDeleteLm( 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        newmonom = pCopy( m );
481        pIncrExp( newmonom, varpermutation[k] );
482        pSetm( newmonom );
483        done= FALSE;
484        while ( list.hasItem() && (done == FALSE) ) {
485            if ( (state= pCmp( list.getItem().monom, newmonom )) < 0 )
486                list++;
487            else done= TRUE;
488        }
489        if ( done == FALSE ) {
490            nlist.append( fglmSelem( newmonom, varpermutation[k] ) );
491            break;
492        }
493        if ( state == 0 ) {
494            list.getItem().newDivisor( varpermutation[k] );
495            pDeleteLm(&newmonom);
496        }
497        else {
498            list.insert( fglmSelem( newmonom, varpermutation[k] ) );
499        }
500        k--;
501    }
502    while ( --k >= 1 ) {
503        newmonom= pCopy( m ); // HIER
504        pIncrExp( newmonom, varpermutation[k] );
505        pSetm( newmonom );
506        nlist.append( fglmSelem( newmonom, varpermutation[k] ) );
507    }
508}
509
510//     if p == pHead( (theIdeal->m)[k] ) return k, 0 otherwise
511//     (Assumes that pLmEqual just checks the leading monomials without
512//      coefficients.)
513int
514fglmSdata::getEdgeNumber( const poly m ) const
515{
516    for ( int k = idelems; k > 0; k-- )
517        if ( pLmEqual( m, (theIdeal->m)[k-1] ) )
518            return k;
519    return 0;
520}
521
522//     Returns the fglmVector v, s.t.
523//        p = v[1]*basis(1) + .. + v[basisSize]*basis(basisSize)
524//     So the size of v depends on the current size of the basis.
525//     Assumes that such an representation exists, i.e. all monoms in p have to be
526//      smaller than basis[basisSize] and that basis[k] < basis[l] for k < l.
527fglmVector
528fglmSdata::getVectorRep( const poly p )
529{
530    fglmVector temp( basisSize );
531    poly m = p;
532    int num = basisSize;
533    while ( m != NULL ) {
534        int comp = pCmp( m, basis[num] );
535        if ( comp == 0 ) {
536            fglmASSERT( num > 0, "Error(1) in fglmSdata::getVectorRep" );
537            number newelem = nCopy( pGetCoeff( m ) );
538            temp.setelem( num, newelem );
539            num--;
540            pIter( m );
541        }
542        else {
543            if ( comp < 0 ) {
544                num--;
545            }
546            else {
547                // This is the place where we can detect if the sourceIdeal
548                // is not reduced. In this case m is not in basis[]. Since basis[]
549                // is ordered this is the case, if and only if basis[i]<m
550                // and basis[j]>m for all j>i
551                _state= FALSE;
552                return temp;
553            }
554        }
555    }
556    return temp;
557}
558
559//     Searches through the border for a monomoial bm which devides m and returns
560//      its normalform in vector representation.
561//     var contains the number of the variable v, s.t. bm = m * v
562fglmVector
563fglmSdata::getBorderDiv( const poly m, int & var ) const
564{
565//     int num2 = borderSize;
566//     while ( num2 > 0 ) {
567//         poly temp = border[num2].monom;
568//         if ( pDivisibleBy( temp, m ) ) {
569//             poly divisor = pDivideM( m, temp );
570//             int var = pIsPurePower( divisor );
571//             if ( (var != 0) && (pGetCoeff( divisor, var ) == 1) ) {
572//                 Print( "poly %s divides poly %s", pString( temp ), pString( m ) );
573//             }
574//         }
575//         num2--;
576//     }
577    int num = borderSize;
578    while ( num > 0 ) {
579        poly temp = border[num].monom;
580        if ( pDivisibleBy( temp, m ) ) {
581            var = pVariables;
582            while ( var > 0 ) {
583                if ( (pGetExp( m, var ) - pGetExp( temp, var )) == 1 )
584                    return border[num].nf;
585                var--;
586            }
587        }
588        num--;
589    }
590    return fglmVector();
591}
592
593void
594internalCalculateFunctionals( const ideal & theIdeal, idealFunctionals & l,
595                              fglmSdata & data )
596{
597
598    // insert pOne() into basis and update the workingList:
599    poly one = pOne();
600    data.newBasisElem( one );
601    data.updateCandidates();
602
603    STICKYPROT(".");
604    while ( data.candidatesLeft() == TRUE ) {
605        fglmSelem candidate = data.nextCandidate();
606        if ( candidate.isBasisOrEdge() == TRUE ) {
607            int edge = data.getEdgeNumber( candidate.monom );
608            if ( edge != 0 ) {
609                // now candidate is an edge, i.e. we know its normalform:
610                // NF(p) = - ( tail(p)/LC(p) )
611                poly nf = data.getSpanPoly( edge );
612                pNorm( nf );
613                pDeleteLm(&nf);  //. deletes the leadingmonomial
614                nf= pNeg( nf );
615                fglmVector nfv = data.getVectorRep( nf );
616                l.insertCols( candidate.divisors, nfv );
617                data.newBorderElem( candidate.monom, nfv );
618                pDelete( &nf );
619                STICKYPROT( "+" );
620            }
621            else {
622                int basis= data.newBasisElem( candidate.monom );
623                data.updateCandidates();
624                l.insertCols( candidate.divisors, basis );
625                STICKYPROT( "." );
626            }
627        }
628        else {
629            int var = 0;
630            fglmVector temp = data.getBorderDiv( candidate.monom, var );
631            fglmASSERT( var > 0, "this should never happen" );
632            fglmVector nfv = l.addCols( var, data.getBasisSize(), temp );
633            data.newBorderElem( candidate.monom, nfv );
634            l.insertCols( candidate.divisors, nfv );
635            STICKYPROT( "-" );
636        }
637        candidate.cleanup();
638    } //. while ( data.candidatesLeft() == TRUE )
639    l.endofConstruction();
640    STICKYPROT2( "\nvdim= %i\n", data.getBasisSize() );
641    return;
642}
643
644//     Calculates the defining Functionals for the ideal "theIdeal" and
645//     returns them in "l".
646//     The ideal has to be zero-dimensional and reduced and has to be a
647//     real subset of the polynomal ring.
648//     In any case it has to be zero-dimensional and minimal (check this
649//      via fglmIdealcheck). Any minimal but not reduced ideal is detected.
650//      In this case it returns FglmNotReduced.
651//     If the base domain is Q, the leading coefficients of the polys
652//     have to be in Z.
653//     returns TRUE if the result is valid, FALSE if theIdeal
654//      is not reduced.
655static BOOLEAN
656CalculateFunctionals( const ideal & theIdeal, idealFunctionals & l )
657{
658    fglmSdata data( theIdeal );
659    internalCalculateFunctionals( theIdeal, l, data );
660    return ( data.state() );
661}
662
663static BOOLEAN
664CalculateFunctionals( const ideal & theIdeal, idealFunctionals & l,
665                      poly & p, fglmVector & v )
666{
667    fglmSdata data( theIdeal );
668    internalCalculateFunctionals( theIdeal, l, data );
669    //    STICKYPROT("Calculating vector rep\n");
670    v = data.getVectorRep( p );
671    // if ( v.isZero() )
672    //   STICKYPROT("vectorrep is 0\n");
673    return ( data.state() );
674}
675
676// ============================================================
677//!      The new basis
678// ============================================================
679
680//     The declaration of class fglmDelem is found in fglm.h
681
682fglmDelem::fglmDelem( poly & m, fglmVector mv, int v ) : v( mv ), insertions( 0 ), var( v )
683{
684    monom= m;
685    m= NULL;
686    for ( int k = pVariables; k > 0; k-- )
687        if ( pGetExp( monom, k ) > 0 )
688            insertions++;
689    // Wir gehen davon aus, dass ein fglmDelem direkt bei der Erzeugung
690    // auch in eine Liste eingefuegt wird. Daher wird hier automatisch
691    // newDivisor aufgerufen ( v teilt ja m )
692    newDivisor();
693}
694
695void
696fglmDelem::cleanup()
697{
698    if ( monom != NULL ) {
699        pDeleteLm(&monom);
700    }
701}
702
703class oldGaussElem
704{
705public:
706    fglmVector v;
707    fglmVector p;
708    number pdenom;
709    number fac;
710
711#ifndef HAVE_EXPLICIT_CONSTR
712    oldGaussElem() : v(), p(), pdenom( NULL ), fac( NULL ) {}
713#endif
714    oldGaussElem( const fglmVector newv, const fglmVector newp, number & newpdenom, number & newfac ) : v( newv ), p( newp ), pdenom( newpdenom ), fac( newfac )
715    {
716        newpdenom= NULL;
717        newfac= NULL;
718    }
719    ~oldGaussElem();
720#ifndef HAVE_EXPLICIT_CONSTR
721    void insertElem( const fglmVector newv, const fglmVector newp, number & newpdenom, number & newfac )
722    {
723        v= newv;
724        p= newp;
725        pdenom= newpdenom;
726        fac= newfac;
727        newpdenom= NULL;
728        newfac= NULL;
729    }
730#endif
731};
732
733oldGaussElem::~oldGaussElem()
734{
735    nDelete( & fac );
736    nDelete( & pdenom );
737}
738
739
740class fglmDdata
741{
742private:
743    int dimen;
744    oldGaussElem * gauss;
745    BOOLEAN * isPivot;  // [1]..[dimen]
746    int * perm;  // [1]..[dimen]
747    int basisSize;  //. the CURRENT basisSize, i.e. basisSize <= dimen
748    polyset basis;  // [1]..[dimen]. The monoms of the new Vectorspace-basis
749
750    int* varpermutation;
751
752    int groebnerBS;
753    int groebnerSize;
754    ideal destId;
755
756    List<fglmDelem> nlist;
757public:
758    fglmDdata( int dimension );
759    ~fglmDdata();
760
761    int getBasisSize() const { return basisSize; }
762    BOOLEAN candidatesLeft() const { return ( nlist.isEmpty() ? FALSE : TRUE ); }
763    fglmDelem nextCandidate();
764    void newBasisElem( poly & m, fglmVector v, fglmVector p, number & denom );
765    void updateCandidates( poly m, const fglmVector v );
766    void newGroebnerPoly( fglmVector & v, poly & p );
767    void gaussreduce( fglmVector & v, fglmVector & p, number & denom );
768    ideal buildIdeal()
769    {
770        idSkipZeroes( destId );
771        return destId;
772    }
773};
774
775fglmDdata::fglmDdata( int dimension )
776{
777    int k;
778    dimen= dimension;
779    basisSize= 0;
780    //. All arrays run from [1]..[dimen], thus omAlloc( dimen + 1 )!
781#ifndef HAVE_EXPLICIT_CONSTR
782    gauss= new oldGaussElem[ dimen+1 ];
783#else
784    gauss= (oldGaussElem *)omAlloc( (dimen+1)*sizeof( oldGaussElem ) );
785#endif
786    isPivot= (BOOLEAN *)omAlloc( (dimen+1)*sizeof( BOOLEAN ) );
787    for ( k= dimen; k > 0; k-- ) isPivot[k]= FALSE;
788    perm= (int *)omAlloc( (dimen+1)*sizeof( int ) );
789    basis= (polyset)omAlloc( (dimen+1)*sizeof( poly ) );
790    varpermutation = (int*)omAlloc( (pVariables+1)*sizeof(int) );
791    // Sort ring variables by increasing values (because of weighted orderings)
792    ideal perm = idMaxIdeal(1);
793    intvec *iv = idSort(perm,TRUE);
794    idDelete(&perm);
795    for(int i = pVariables; i > 0; i--) varpermutation[pVariables+1-i] = (*iv)[i-1];
796    delete iv;
797
798    groebnerBS= 16;
799    groebnerSize= 0;
800    destId= idInit( groebnerBS, 1 );
801}
802
803fglmDdata::~fglmDdata()
804{
805  // STICKYPROT2("dimen= %i", dimen);
806  // STICKYPROT2("basisSize= %i", basisSize);
807  //    fglmASSERT( dimen == basisSize, "Es wurden nicht alle BasisElemente gefunden!" );
808    int k;
809#ifndef HAVE_EXPLICIT_CONSTR
810    delete [] gauss;
811#else
812    // use basisSize instead of dimen because of fglmquot!
813    for ( k= basisSize; k > 0; k-- )
814        gauss[k].~oldGaussElem();
815    omFreeSize( (ADDRESS)gauss, (dimen+1)*sizeof( oldGaussElem ) );
816#endif
817    omFreeSize( (ADDRESS)isPivot, (dimen+1)*sizeof( BOOLEAN ) );
818    omFreeSize( (ADDRESS)perm, (dimen+1)*sizeof( int ) );
819    // use basisSize instead of dimen because of fglmquot!
820    //. Remember: There is no poly in basis[0], thus k > 0
821    for ( k= basisSize; k > 0; k-- )
822        pLmDelete( basis[k]);
823    omFreeSize( (ADDRESS)basis, (dimen+1)*sizeof( poly ) );
824    omFreeSize( (ADDRESS)varpermutation, (pVariables+1)*sizeof(int) );
825}
826
827fglmDelem
828fglmDdata::nextCandidate()
829{
830    fglmDelem result = nlist.getFirst();
831    nlist.removeFirst();
832    return result;
833}
834
835void
836fglmDdata::newBasisElem( poly & m, fglmVector v, fglmVector p, number & denom )
837{
838// inserts m as a new basis monom. m is NOT copied but directly inserted.
839// returns m=NULL to indicate, that now basis is oweing m.
840    basisSize++;
841    basis[basisSize]= m;
842    m= NULL;
843    int k= 1;
844    while ( nIsZero(v.getconstelem(k)) || isPivot[k] ) {
845        k++;
846    }
847    fglmASSERT( k <= dimen, "Error(1) in fglmDdata::pivot-search");
848    number pivot= v.getconstelem( k );
849    int pivotcol = k;
850    k++;
851    while ( k <= dimen ) {
852        if ( ! nIsZero( v.getconstelem(k) ) && ! isPivot[k] ) {
853            if ( nGreater( v.getconstelem( k ), pivot ) ) {
854                pivot= v.getconstelem( k );
855                pivotcol= k;
856            }
857        }
858        k++;
859    }
860    fglmASSERT( ! nIsZero( pivot ), "Error(2) fglmDdata::Pivotelement ist Null" );
861    isPivot[ pivotcol ]= TRUE;
862    perm[basisSize]= pivotcol;
863
864    pivot= nCopy( v.getconstelem( pivotcol ) );
865#ifndef HAVE_EXPLICIT_CONSTR
866    gauss[basisSize].insertElem( v, p, denom, pivot );
867#else
868    gauss[basisSize].oldGaussElem( v, p, denom, pivot );
869#endif
870}
871
872void
873fglmDdata::updateCandidates( poly m, const fglmVector v )
874{
875    ListIterator<fglmDelem> list = nlist;
876    poly newmonom = NULL;
877    int k = pVariables;
878    BOOLEAN done = FALSE;
879    int state = 0;
880    while ( k >= 1 ) {
881        newmonom = pCopy( m );
882        pIncrExp( newmonom, varpermutation[k] );
883        pSetm( newmonom );
884        done= FALSE;
885        while ( list.hasItem() && (done == FALSE) ) {
886            if ( (state= pCmp( list.getItem().monom, newmonom )) < 0 )
887                list++;
888            else done= TRUE;
889        }
890        if ( done == FALSE ) {
891            nlist.append( fglmDelem( newmonom, v, k ) );
892            break;
893        }
894        if ( state == 0 ) {
895            list.getItem().newDivisor();
896            pDeleteLm( & newmonom );
897        }
898        else {
899            list.insert( fglmDelem( newmonom, v, k ) );
900        }
901        k--;
902    }
903    while ( --k >= 1 ) {
904        newmonom= pCopy( m );
905        pIncrExp( newmonom, varpermutation[k] );
906        pSetm( newmonom );
907        nlist.append( fglmDelem( newmonom, v, k ) );
908    }
909}
910
911void
912fglmDdata::newGroebnerPoly( fglmVector & p, poly & m )
913// Inserts gp = p[1]*basis(1)+..+p[basisSize]*basis(basisSize)+p[basisSize+1]*m as
914//  a new groebner polynomial for the ideal.
915// All elements (monomials and coefficients) of gp are copied, instead of m.
916// Assumes that p.length() == basisSize+1.
917{
918    //. Baue das Polynom von oben nach unten:
919    fglmASSERT( p.size() == basisSize+1, "GP::newGroebnerPoly: p has wrong size" );
920    int k;
921    poly result = m;
922    poly temp = result;
923    m= NULL;
924    if ( nGetChar() > 0 ) {
925        number lead = nCopy( p.getconstelem( basisSize+1 ) );
926        p /= lead;
927        nDelete( & lead );
928    }
929    if ( nGetChar() == 0 ) {
930        number gcd= p.gcd();
931        fglmASSERT( ! nIsZero( gcd ), "FATAL: gcd and thus p is zero" );
932        if ( ! nIsOne( gcd ) )
933            p /= gcd;
934        nDelete( & gcd );
935    }
936    pSetCoeff( result, nCopy( p.getconstelem( basisSize+1 ) ) );
937    for ( k= basisSize; k > 0; k-- ) {
938        if ( ! nIsZero( p.getconstelem( k ) ) ) {
939            temp->next= pCopy( basis[k] );
940            pIter( temp );
941            pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
942        }
943    }
944    pSetm( result );
945    if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
946    if ( groebnerSize == IDELEMS( destId ) ) {
947        pEnlargeSet( & destId->m, IDELEMS( destId ), groebnerBS );
948        IDELEMS( destId )+= groebnerBS;
949    }
950    (destId->m)[groebnerSize]= result;
951    groebnerSize++;
952}
953
954void
955fglmDdata::gaussreduce( fglmVector & v, fglmVector & p, number & pdenom )
956{
957    int k;
958    number fac1, fac2;
959    number temp;
960    fglmASSERT( pdenom == NULL, "pdenom in gaussreduce should be NULL" );
961    pdenom= nInit( 1 );
962    number vdenom = v.clearDenom();
963    if ( ! nIsOne( vdenom ) && ! nIsZero( vdenom ) ) {
964        p.setelem( p.size(), vdenom );
965    }
966    else {
967        nDelete( &vdenom );
968    }
969    number gcd = v.gcd();
970    if ( ! nIsOne( gcd ) && ! nIsZero( gcd ) ) {
971        v /= gcd;
972        number temp= nMult( pdenom, gcd );
973        nDelete( &pdenom );
974        pdenom= temp;
975    }
976    nDelete( & gcd );
977
978    for ( k= 1; k <= basisSize; k++ ) {
979
980        if ( ! v.elemIsZero( perm[k] ) ) {
981            fac1= gauss[k].fac;
982            fac2= nCopy( v.getconstelem( perm[k] ) );
983            v.nihilate( fac1, fac2, gauss[k].v );
984            fac1= nMult( fac1, gauss[k].pdenom );
985            temp= nMult( fac2, pdenom );
986            nDelete( &fac2 );
987            fac2= temp;
988              p.nihilate( fac1, fac2, gauss[k].p );
989            temp= nMult( pdenom, gauss[k].pdenom );
990            nDelete( &pdenom );
991            pdenom= temp;
992
993            nDelete( & fac1 );
994            nDelete( & fac2 );
995            number gcd = v.gcd();
996            if ( ! nIsOne( gcd ) && ! nIsZero( gcd ) ) {
997                v /= gcd;
998                number temp= nMult( pdenom, gcd );
999                nDelete( &pdenom );
1000                pdenom= temp;
1001            }
1002            nDelete( & gcd );
1003            gcd= p.gcd();
1004            temp= nGcd( pdenom, gcd, currRing );
1005            nDelete( &gcd );
1006            gcd= temp;
1007            if ( ! nIsZero( gcd ) && ! nIsOne( gcd ) ) {
1008                p /= gcd;
1009                temp= nDiv( pdenom, gcd );
1010                nDelete( & pdenom );
1011                pdenom= temp;
1012                nNormalize( pdenom );
1013            }
1014            nDelete( & gcd );
1015        }
1016    }
1017}
1018
1019static ideal
1020GroebnerViaFunctionals( const idealFunctionals & l,
1021                        fglmVector iv = fglmVector() )
1022// If iv is zero, calculates the groebnerBasis for the ideal which is
1023// defined by l.
1024// If iv is not zero, then the groebnerBasis if i:p is calculated where
1025// i is defined by l and iv is the vector-representation of nf(p) wrt. i
1026// The dimension of l has to be finite.
1027// The result is in reduced form.
1028{
1029    fglmDdata data( l.dimen() );
1030
1031    // insert pOne() and update workinglist according to iv:
1032    fglmVector initv;
1033    if ( iv.isZero() ) {
1034      // STICKYPROT("initv is zero\n");
1035      initv = fglmVector( l.dimen(), 1 );
1036    }
1037    else {
1038      // STICKYPROT("initv is not zero\n");
1039      initv = iv;
1040    }
1041     
1042    poly one = pOne();
1043    data.updateCandidates( one, initv );
1044    number nOne = nInit( 1 );
1045    data.newBasisElem( one, initv, fglmVector( 1, 1 ), nOne );
1046    STICKYPROT( "." );
1047    while ( data.candidatesLeft() == TRUE ) {
1048        fglmDelem candidate = data.nextCandidate();
1049        if ( candidate.isBasisOrEdge() == TRUE ) {
1050            // Now we have the chance to find a new groebner polynomial
1051
1052            // v is the vector-representation of candidate.monom
1053            // some elements of v are zeroed in data.gaussreduce(). Which
1054            // ones and how this was done is stored in p.
1055            // originalV containes the unchanged v, which is later inserted
1056            // into the working list (via data.updateCandidates().
1057            fglmVector v = l.multiply( candidate.v, candidate.var );
1058            fglmVector originalV = v;
1059            fglmVector p( data.getBasisSize()+1, data.getBasisSize()+1 );
1060            number pdenom = NULL;
1061            data.gaussreduce( v, p, pdenom );
1062            if ( v.isZero() ) {
1063                // Now v is linear dependend to the already found basis elements.
1064                // This means that v (rsp. candidate.monom) is the leading
1065                // monomial of the next groebner-basis polynomial.
1066                data.newGroebnerPoly( p, candidate.monom );
1067                nDelete( & pdenom );
1068                STICKYPROT( "+" );
1069            }
1070            else {
1071                // no linear dependence could be found, so v ( rsp. monom )
1072                // is a basis monomial. We store the zeroed version ( i.e. v
1073                // and not originalV ) as well as p, the denomiator and all
1074                // the other stuff.
1075                // erst updateCandidates, dann newBasisELem!!!
1076                data.updateCandidates( candidate.monom, originalV );
1077                data.newBasisElem( candidate.monom, v, p, pdenom );
1078                STICKYPROT( "." );
1079            }
1080        }
1081        else {
1082            STICKYPROT( "-" );
1083            candidate.cleanup();
1084        }
1085    }  //. while data.candidatesLeft()
1086    STICKYPROT( "\n" );
1087    return ( data.buildIdeal() );
1088}
1089//<-
1090
1091static ideal
1092FindUnivariatePolys( const idealFunctionals & l )
1093{
1094    fglmVector v;
1095    fglmVector p;
1096    ideal destIdeal = idInit( pVariables, 1 );
1097
1098    int i;
1099    BOOLEAN isZero;
1100    int *varpermutation = (int*)omAlloc( (pVariables+1)*sizeof(int) );
1101    ideal perm = idMaxIdeal(1);
1102    intvec *iv = idSort(perm,TRUE);
1103    idDelete(&perm);
1104    for(i = pVariables; i > 0; i--) varpermutation[pVariables+1-i] = (*iv)[i-1];
1105    delete iv;
1106     
1107    for (i= 1; i <= pVariables; i++ )
1108    {
1109        // main loop
1110        STICKYPROT2( "(%i)", i /*varpermutation[i]*/);
1111        gaussReducer gauss( l.dimen() );
1112        isZero= FALSE;
1113        v= fglmVector( l.dimen(), 1 );
1114        while ( isZero == FALSE ) {
1115            if ( (isZero= gauss.reduce( v )) == TRUE ) {
1116                STICKYPROT( "+" );
1117                p= gauss.getDependence();
1118                number gcd= p.gcd();
1119                if ( ! nIsOne( gcd ) ) {
1120                    p /= gcd;
1121                }
1122                nDelete( & gcd );
1123                int k;
1124                poly temp = NULL;
1125                poly result;
1126                for ( k= p.size(); k > 0; k-- ) {
1127                    number n = nCopy( p.getconstelem( k ) );
1128                    if ( ! nIsZero( n ) ) {
1129                        if ( temp == NULL ) {
1130                            result= pOne();
1131                            temp= result;
1132                        }
1133                        else {
1134                            temp->next= pOne();
1135                            pIter( temp );
1136                        }
1137                        pSetCoeff( temp, n );
1138                        pSetExp( temp, i /*varpermutation[i]*/, k-1 );
1139                        pSetm( temp );
1140                    }
1141                }
1142                if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
1143                (destIdeal->m)[i-1]= result;
1144            }
1145            else {
1146                STICKYPROT( "." );
1147                gauss.store();
1148                v= l.multiply( v, i /*varpermutation[i]*/ );
1149            }
1150        }
1151    }
1152    STICKYPROT( "\n" );
1153    omFreeSize( (ADDRESS)varpermutation, (pVariables+1)*sizeof(int) );
1154    return destIdeal;
1155}
1156
1157// for a descritption of the parameters see fglm.h
1158BOOLEAN
1159fglmzero( idhdl sourceRingHdl, ideal & sourceIdeal, idhdl destRingHdl, ideal & destIdeal, BOOLEAN switchBack, BOOLEAN deleteIdeal )
1160{
1161    idhdl initialRingHdl = currRingHdl;
1162    BOOLEAN fglmok;
1163
1164    if ( currRingHdl != sourceRingHdl )
1165        rSetHdl( sourceRingHdl );
1166    idealFunctionals L( 100, pVariables );
1167    fglmok = CalculateFunctionals( sourceIdeal, L );
1168    if ( deleteIdeal == TRUE )
1169        idDelete( & sourceIdeal );
1170    rSetHdl( destRingHdl );
1171    if ( fglmok == TRUE ) {
1172        L.map( IDRING( sourceRingHdl ) );
1173        destIdeal= GroebnerViaFunctionals( L );
1174    }
1175    if ( (switchBack == TRUE) && (currRingHdl != initialRingHdl) )
1176        rSetHdl( initialRingHdl );
1177    return fglmok;
1178}
1179
1180BOOLEAN
1181fglmquot( ideal sourceIdeal, poly quot, ideal & destIdeal)
1182{
1183    BOOLEAN fglmok;
1184    fglmVector v;
1185
1186    idealFunctionals L( 100, pVariables );
1187    // STICKYPROT("calculating normal form\n");
1188    // poly p = kNF( sourceIdeal, currRing->qideal, quot );
1189    // STICKYPROT("calculating functionals\n");
1190    fglmok = CalculateFunctionals( sourceIdeal, L, quot, v );
1191    if ( fglmok == TRUE ) {
1192      // STICKYPROT("calculating groebner basis\n");
1193        destIdeal= GroebnerViaFunctionals( L, v );
1194    }
1195    return fglmok;
1196}
1197
1198BOOLEAN
1199FindUnivariateWrapper( ideal source, ideal & destIdeal )
1200{
1201    BOOLEAN fglmok;
1202
1203    idealFunctionals L( 100, pVariables );
1204    fglmok = CalculateFunctionals( source, L );
1205    if ( fglmok == TRUE ) {
1206        destIdeal= FindUnivariatePolys( L );
1207        return TRUE;
1208    }
1209    else
1210        return FALSE;
1211}
1212
1213
1214#endif
1215
1216// ----------------------------------------------------------------------------
1217// Local Variables: ***
1218// compile-command: "make Singular" ***
1219// page-delimiter: "^\\(\\|//!\\)" ***
1220// fold-internal-margins: nil ***
1221// End: ***
Note: See TracBrowser for help on using the repository browser.