source: git/Singular/fglmzero.cc @ c232af

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