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