[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 number-vectors for the fglm algorithm. |
---|
| 10 | * (See fglm.cc). Based on a letter-envelope implementation, mainly |
---|
| 11 | * written to be used by the fglm algorithm. Hence they are |
---|
| 12 | * specialized for this purpose. |
---|
| 13 | */ |
---|
| 14 | |
---|
[599326] | 15 | #include <kernel/mod2.h> |
---|
[35aab3] | 16 | |
---|
[6eccc9] | 17 | #ifdef HAVE_FACTORY |
---|
[b1dfaf] | 18 | #include <omalloc/omalloc.h> |
---|
[599326] | 19 | #include <kernel/structs.h> |
---|
[0f401f] | 20 | #include <coeffs/numbers.h> |
---|
[599326] | 21 | #include <kernel/fglm.h> |
---|
| 22 | #include <kernel/fglmvec.h> |
---|
[35aab3] | 23 | |
---|
| 24 | #define PROT(msg) |
---|
| 25 | #define STICKYPROT(msg) if (BTEST1(OPT_PROT)) Print(msg) |
---|
| 26 | #define PROT2(msg,arg) |
---|
| 27 | #define STICKYPROT2(msg,arg) if (BTEST1(OPT_PROT)) Print(msg,arg) |
---|
| 28 | #define fglmASSERT(ignore1,ignore2) |
---|
| 29 | |
---|
| 30 | class fglmVectorRep |
---|
| 31 | { |
---|
| 32 | private: |
---|
| 33 | int ref_count; |
---|
| 34 | int N; |
---|
| 35 | number * elems; |
---|
| 36 | public: |
---|
| 37 | fglmVectorRep() : ref_count(1), N(0), elems(0) {} |
---|
| 38 | fglmVectorRep( int n, number * e ) : ref_count(1), N(n), elems(e) {} |
---|
| 39 | fglmVectorRep( int n ) : ref_count(1), N(n) |
---|
| 40 | { |
---|
| 41 | fglmASSERT( N >= 0, "illegal Vector representation" ); |
---|
| 42 | if ( N == 0 ) |
---|
| 43 | elems= 0; |
---|
| 44 | else { |
---|
| 45 | elems= (number *)omAlloc( N*sizeof( number ) ); |
---|
| 46 | for ( int i= N-1; i >= 0; i-- ) |
---|
| 47 | elems[i]= nInit( 0 ); |
---|
| 48 | } |
---|
| 49 | } |
---|
| 50 | ~fglmVectorRep() |
---|
| 51 | { |
---|
| 52 | if ( N > 0 ) { |
---|
| 53 | for ( int i= N-1; i >= 0; i-- ) |
---|
| 54 | nDelete( elems + i ); |
---|
| 55 | omFreeSize( (ADDRESS)elems, N*sizeof( number ) ); |
---|
| 56 | } |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | fglmVectorRep* clone() const |
---|
| 60 | { |
---|
| 61 | if ( N > 0 ) { |
---|
| 62 | number * elems_clone; |
---|
| 63 | elems_clone= (number *)omAlloc( N*sizeof( number ) ); |
---|
| 64 | for ( int i= N-1; i >= 0; i-- ) |
---|
| 65 | elems_clone[i] = nCopy( elems[i] ); |
---|
| 66 | return new fglmVectorRep( N, elems_clone ); |
---|
| 67 | } else |
---|
| 68 | return new fglmVectorRep( N, 0 ); |
---|
| 69 | } |
---|
| 70 | BOOLEAN deleteObject() { return --ref_count == 0; } |
---|
| 71 | fglmVectorRep * copyObject() { ref_count++; return this; } |
---|
| 72 | int refcount() const { return ref_count; } |
---|
| 73 | BOOLEAN isUnique() const { return ref_count == 1; } |
---|
| 74 | |
---|
| 75 | int size() const { return N; } |
---|
| 76 | int isZero() const |
---|
| 77 | { |
---|
| 78 | int k; |
---|
| 79 | for ( k= N; k > 0; k-- ) |
---|
| 80 | if ( ! nIsZero( getconstelem( k ) ) ) |
---|
| 81 | return 0; |
---|
| 82 | return 1; |
---|
| 83 | } |
---|
| 84 | int numNonZeroElems() const |
---|
| 85 | { |
---|
| 86 | int num = 0; |
---|
| 87 | int k; |
---|
| 88 | for ( k= N; k > 0; k-- ) |
---|
| 89 | if ( ! nIsZero( getconstelem( k ) ) ) num++; |
---|
| 90 | return num; |
---|
| 91 | } |
---|
| 92 | void setelem( int i, number n ) |
---|
| 93 | { |
---|
| 94 | fglmASSERT( 0 < i && i <= N, "setelem: wrong index" ); |
---|
| 95 | nDelete( elems + i-1 ); |
---|
| 96 | elems[i-1]= n; |
---|
| 97 | } |
---|
| 98 | number ejectelem( int i, number n ) |
---|
| 99 | { |
---|
| 100 | fglmASSERT( isUnique(), "should only be called if unique!" ); |
---|
| 101 | number temp= elems[i-1]; |
---|
| 102 | elems[i-1]= n; |
---|
| 103 | return temp; |
---|
| 104 | } |
---|
| 105 | number & getelem( int i ) |
---|
| 106 | { |
---|
| 107 | fglmASSERT( 0 < i && i <= N, "getelem: wrong index" ); |
---|
| 108 | return elems[i-1]; |
---|
| 109 | } |
---|
| 110 | const number getconstelem( int i) const |
---|
| 111 | { |
---|
| 112 | fglmASSERT( 0 < i && i <= N, "getconstelem: wrong index" ); |
---|
| 113 | return elems[i-1]; |
---|
| 114 | } |
---|
| 115 | friend class fglmVector; |
---|
| 116 | }; |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | ///-------------------------------------------------------------------------------- |
---|
| 120 | /// Implementation of class fglmVector |
---|
| 121 | ///-------------------------------------------------------------------------------- |
---|
| 122 | |
---|
| 123 | fglmVector::fglmVector( fglmVectorRep * r ) : rep(r) {} |
---|
| 124 | |
---|
| 125 | fglmVector::fglmVector() : rep( new fglmVectorRep() ) {} |
---|
| 126 | |
---|
| 127 | fglmVector::fglmVector( int size ) : rep( new fglmVectorRep( size ) ) {} |
---|
| 128 | |
---|
| 129 | fglmVector::fglmVector( int size, int basis ) : rep( new fglmVectorRep( size ) ) |
---|
| 130 | { |
---|
| 131 | rep->setelem( basis, nInit( 1 ) ); |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | fglmVector::fglmVector( const fglmVector & v ) |
---|
| 135 | { |
---|
| 136 | rep= v.rep->copyObject(); |
---|
| 137 | } |
---|
| 138 | |
---|
| 139 | fglmVector::~fglmVector() |
---|
| 140 | { |
---|
| 141 | if ( rep->deleteObject() ) |
---|
| 142 | delete rep; |
---|
| 143 | } |
---|
| 144 | |
---|
| 145 | #ifndef HAVE_EXPLICIT_CONSTR |
---|
| 146 | void |
---|
| 147 | fglmVector::mac_constr( const fglmVector & v) |
---|
| 148 | { |
---|
| 149 | rep= v.rep->copyObject(); |
---|
| 150 | } |
---|
| 151 | |
---|
| 152 | void |
---|
| 153 | fglmVector::mac_constr_i( int size ) |
---|
| 154 | { |
---|
| 155 | rep= new fglmVectorRep( size ); |
---|
| 156 | } |
---|
| 157 | |
---|
| 158 | void |
---|
| 159 | fglmVector::clearelems() |
---|
| 160 | { |
---|
| 161 | if ( rep->deleteObject() ) |
---|
| 162 | delete rep; |
---|
| 163 | } |
---|
| 164 | #endif |
---|
| 165 | |
---|
| 166 | void |
---|
| 167 | fglmVector::makeUnique() |
---|
| 168 | { |
---|
| 169 | if ( rep->refcount() != 1 ) { |
---|
| 170 | rep->deleteObject(); |
---|
| 171 | rep= rep->clone(); |
---|
| 172 | } |
---|
| 173 | } |
---|
| 174 | |
---|
| 175 | int |
---|
| 176 | fglmVector::size() const |
---|
| 177 | { |
---|
| 178 | return rep->size(); |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | int |
---|
| 182 | fglmVector::numNonZeroElems() const |
---|
| 183 | { |
---|
| 184 | return rep->numNonZeroElems(); |
---|
| 185 | } |
---|
| 186 | |
---|
| 187 | void |
---|
| 188 | fglmVector::nihilate( const number fac1, const number fac2, const fglmVector v ) |
---|
| 189 | { |
---|
| 190 | int i; |
---|
| 191 | int vsize= v.size(); |
---|
| 192 | number term1, term2; |
---|
| 193 | fglmASSERT( vsize <= rep->size(), "v has to be smaller oder equal" ); |
---|
| 194 | if ( rep->isUnique() ) { |
---|
| 195 | for ( i= vsize; i > 0; i-- ) { |
---|
| 196 | term1= nMult( fac1, rep->getconstelem( i ) ); |
---|
| 197 | term2= nMult( fac2, v.rep->getconstelem( i ) ); |
---|
| 198 | rep->setelem( i, nSub( term1, term2 ) ); |
---|
| 199 | nDelete( &term1 ); |
---|
| 200 | nDelete( &term2 ); |
---|
| 201 | } |
---|
| 202 | for ( i= rep->size(); i > vsize; i-- ) { |
---|
| 203 | rep->setelem( i, nMult( fac1, rep->getconstelem( i ) ) ); |
---|
| 204 | } |
---|
| 205 | } |
---|
| 206 | else |
---|
| 207 | { |
---|
| 208 | number* newelems; |
---|
| 209 | newelems= (number *)omAlloc( rep->size()*sizeof( number ) ); |
---|
| 210 | for ( i= vsize; i > 0; i-- ) { |
---|
| 211 | term1= nMult( fac1, rep->getconstelem( i ) ); |
---|
| 212 | term2= nMult( fac2, v.rep->getconstelem( i ) ); |
---|
| 213 | newelems[i-1]= nSub( term1, term2 ); |
---|
| 214 | nDelete( &term1 ); |
---|
| 215 | nDelete( &term2 ); |
---|
| 216 | } |
---|
| 217 | for ( i= rep->size(); i > vsize; i-- ) { |
---|
| 218 | newelems[i-1]= nMult( fac1, rep->getconstelem( i ) ); |
---|
| 219 | } |
---|
| 220 | rep->deleteObject(); |
---|
| 221 | rep= new fglmVectorRep( rep->size(), newelems ); |
---|
| 222 | } |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | fglmVector & |
---|
| 226 | fglmVector::operator = ( const fglmVector & v ) |
---|
| 227 | { |
---|
| 228 | if ( this != &v ) { |
---|
| 229 | if ( rep->deleteObject() ) |
---|
| 230 | delete rep; |
---|
| 231 | rep = v.rep->copyObject(); |
---|
| 232 | } |
---|
| 233 | return *this; |
---|
| 234 | } |
---|
| 235 | |
---|
| 236 | int |
---|
| 237 | fglmVector::operator == ( const fglmVector & v ) |
---|
| 238 | { |
---|
| 239 | if ( rep->size() == v.rep->size() ) { |
---|
| 240 | if ( rep == v.rep ) return 1; |
---|
| 241 | else { |
---|
| 242 | int i; |
---|
| 243 | for ( i= rep->size(); i > 0; i-- ) |
---|
| 244 | if ( ! nEqual( rep->getconstelem( i ), v.rep->getconstelem( i ) ) ) |
---|
| 245 | return 0; |
---|
| 246 | return 1; |
---|
| 247 | } |
---|
| 248 | } |
---|
| 249 | return 0; |
---|
| 250 | } |
---|
| 251 | |
---|
| 252 | int |
---|
| 253 | fglmVector::operator != ( const fglmVector & v ) |
---|
| 254 | { |
---|
| 255 | return !( *this == v ); |
---|
| 256 | } |
---|
| 257 | |
---|
| 258 | int |
---|
| 259 | fglmVector::isZero() |
---|
| 260 | { |
---|
| 261 | return rep->isZero(); |
---|
| 262 | } |
---|
| 263 | |
---|
| 264 | int |
---|
| 265 | fglmVector::elemIsZero( int i ) |
---|
| 266 | { |
---|
| 267 | return nIsZero( rep->getconstelem( i ) ); |
---|
| 268 | } |
---|
| 269 | |
---|
| 270 | fglmVector & |
---|
| 271 | fglmVector::operator += ( const fglmVector & v ) |
---|
| 272 | { |
---|
| 273 | fglmASSERT( size() == v.size(), "incompatible vectors" ); |
---|
| 274 | // ACHTUNG : Das Verhalten hier mit gcd genau ueberpruefen! |
---|
| 275 | int i; |
---|
| 276 | if ( rep->isUnique() ) { |
---|
| 277 | for ( i= rep->size(); i > 0; i-- ) |
---|
| 278 | rep->setelem( i, nAdd( rep->getconstelem( i ), v.rep->getconstelem( i ) ) ); |
---|
| 279 | } |
---|
| 280 | else |
---|
| 281 | { |
---|
| 282 | int n = rep->size(); |
---|
| 283 | number* newelems; |
---|
| 284 | newelems= (number *)omAlloc( n*sizeof( number ) ); |
---|
| 285 | for ( i= n; i > 0; i-- ) |
---|
| 286 | newelems[i-1]= nAdd( rep->getconstelem( i ), v.rep->getconstelem( i ) ); |
---|
| 287 | rep->deleteObject(); |
---|
| 288 | rep= new fglmVectorRep( n, newelems ); |
---|
| 289 | } |
---|
| 290 | return *this; |
---|
| 291 | } |
---|
| 292 | |
---|
| 293 | fglmVector & |
---|
| 294 | fglmVector::operator -= ( const fglmVector & v ) |
---|
| 295 | { |
---|
| 296 | fglmASSERT( size() == v.size(), "incompatible vectors" ); |
---|
| 297 | int i; |
---|
| 298 | if ( rep->isUnique() ) { |
---|
| 299 | for ( i= rep->size(); i > 0; i-- ) |
---|
| 300 | rep->setelem( i, nSub( rep->getconstelem( i ), v.rep->getconstelem( i ) ) ); |
---|
| 301 | } |
---|
| 302 | else |
---|
| 303 | { |
---|
| 304 | int n = rep->size(); |
---|
| 305 | number* newelems; |
---|
| 306 | newelems= (number *)omAlloc( n*sizeof( number ) ); |
---|
| 307 | for ( i= n; i > 0; i-- ) |
---|
| 308 | newelems[i-1]= nSub( rep->getconstelem( i ), v.rep->getconstelem( i ) ); |
---|
| 309 | rep->deleteObject(); |
---|
| 310 | rep= new fglmVectorRep( n, newelems ); |
---|
| 311 | } |
---|
| 312 | return *this; |
---|
| 313 | } |
---|
| 314 | |
---|
| 315 | fglmVector & |
---|
| 316 | fglmVector::operator *= ( const number & n ) |
---|
| 317 | { |
---|
| 318 | int s= rep->size(); |
---|
| 319 | int i; |
---|
| 320 | if ( ! rep->isUnique() ) { |
---|
| 321 | number * temp; |
---|
| 322 | temp= (number *)omAlloc( s*sizeof( number ) ); |
---|
| 323 | for ( i= s; i > 0; i-- ) |
---|
| 324 | temp[i-1]= nMult( rep->getconstelem( i ), n ); |
---|
| 325 | rep->deleteObject(); |
---|
| 326 | rep= new fglmVectorRep( s, temp ); |
---|
| 327 | } |
---|
| 328 | else |
---|
| 329 | { |
---|
| 330 | for (i= s; i > 0; i-- ) |
---|
| 331 | rep->setelem( i, nMult( rep->getconstelem( i ), n ) ); |
---|
| 332 | } |
---|
| 333 | return *this; |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | fglmVector & |
---|
| 337 | fglmVector::operator /= ( const number & n ) |
---|
| 338 | { |
---|
| 339 | int s= rep->size(); |
---|
| 340 | int i; |
---|
| 341 | if ( ! rep->isUnique() ) { |
---|
| 342 | number * temp; |
---|
| 343 | temp= (number *)omAlloc( s*sizeof( number ) ); |
---|
| 344 | for ( i= s; i > 0; i-- ) { |
---|
| 345 | temp[i-1]= nDiv( rep->getconstelem( i ), n ); |
---|
| 346 | nNormalize( temp[i-1] ); |
---|
| 347 | } |
---|
| 348 | rep->deleteObject(); |
---|
| 349 | rep= new fglmVectorRep( s, temp ); |
---|
| 350 | } |
---|
| 351 | else |
---|
| 352 | { |
---|
| 353 | for (i= s; i > 0; i-- ) { |
---|
| 354 | rep->setelem( i, nDiv( rep->getconstelem( i ), n ) ); |
---|
| 355 | nNormalize( rep->getelem( i ) ); |
---|
| 356 | } |
---|
| 357 | } |
---|
| 358 | return *this; |
---|
| 359 | } |
---|
| 360 | |
---|
| 361 | fglmVector |
---|
| 362 | operator - ( const fglmVector & v ) |
---|
| 363 | { |
---|
| 364 | fglmVector temp( v.size() ); |
---|
| 365 | int i; |
---|
| 366 | number n ; |
---|
| 367 | for ( i= v.size(); i > 0; i-- ) { |
---|
| 368 | n= nCopy( v.getconstelem( i ) ); |
---|
| 369 | n= nNeg( n ); |
---|
| 370 | temp.setelem( i, n ); |
---|
| 371 | } |
---|
| 372 | return temp; |
---|
| 373 | } |
---|
| 374 | |
---|
| 375 | fglmVector |
---|
| 376 | operator + ( const fglmVector & lhs, const fglmVector & rhs ) |
---|
| 377 | { |
---|
| 378 | fglmVector temp= lhs; |
---|
| 379 | temp+= rhs; |
---|
| 380 | return temp; |
---|
| 381 | } |
---|
| 382 | |
---|
| 383 | fglmVector |
---|
| 384 | operator - ( const fglmVector & lhs, const fglmVector & rhs ) |
---|
| 385 | { |
---|
| 386 | fglmVector temp= lhs; |
---|
| 387 | temp-= rhs; |
---|
| 388 | return temp; |
---|
| 389 | } |
---|
| 390 | |
---|
| 391 | fglmVector |
---|
| 392 | operator * ( const fglmVector & v, const number n ) |
---|
| 393 | { |
---|
| 394 | fglmVector temp= v; |
---|
| 395 | temp*= n; |
---|
| 396 | return temp; |
---|
| 397 | } |
---|
| 398 | |
---|
| 399 | fglmVector |
---|
| 400 | operator * ( const number n, const fglmVector & v ) |
---|
| 401 | { |
---|
| 402 | fglmVector temp= v; |
---|
| 403 | temp*= n; |
---|
| 404 | return temp; |
---|
| 405 | } |
---|
| 406 | |
---|
| 407 | number & |
---|
| 408 | fglmVector::getelem( int i ) |
---|
| 409 | { |
---|
| 410 | makeUnique(); |
---|
| 411 | return rep->getelem( i ); |
---|
| 412 | } |
---|
| 413 | |
---|
[16fa35] | 414 | const number |
---|
[35aab3] | 415 | fglmVector::getconstelem( int i ) const |
---|
| 416 | { |
---|
| 417 | return rep->getconstelem( i ); |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | void |
---|
| 421 | fglmVector::setelem( int i, number & n ) |
---|
| 422 | { |
---|
| 423 | makeUnique(); |
---|
| 424 | rep->setelem( i, n ); |
---|
| 425 | n= nNULL; |
---|
| 426 | } |
---|
| 427 | |
---|
| 428 | number |
---|
| 429 | fglmVector::gcd() const |
---|
| 430 | { |
---|
| 431 | int i= rep->size(); |
---|
| 432 | BOOLEAN found = FALSE; |
---|
| 433 | BOOLEAN gcdIsOne = FALSE; |
---|
| 434 | number theGcd; |
---|
| 435 | number current; |
---|
| 436 | while( i > 0 && ! found ) { |
---|
| 437 | current= rep->getconstelem( i ); |
---|
| 438 | if ( ! nIsZero( current ) ) { |
---|
| 439 | theGcd= nCopy( current ); |
---|
| 440 | found= TRUE; |
---|
| 441 | if ( ! nGreaterZero( theGcd ) ) { |
---|
| 442 | theGcd= nNeg( theGcd ); |
---|
| 443 | } |
---|
| 444 | if ( nIsOne( theGcd ) ) gcdIsOne= TRUE; |
---|
| 445 | } |
---|
| 446 | i--; |
---|
| 447 | } |
---|
| 448 | if ( found ) { |
---|
| 449 | while ( i > 0 && ! gcdIsOne ) { |
---|
| 450 | current= rep->getconstelem( i ); |
---|
| 451 | if ( ! nIsZero( current ) ) { |
---|
| 452 | number temp= nGcd( theGcd, current, currRing ); |
---|
| 453 | nDelete( &theGcd ); |
---|
| 454 | theGcd= temp; |
---|
| 455 | if ( nIsOne( theGcd ) ) gcdIsOne= TRUE; |
---|
| 456 | } |
---|
| 457 | i--; |
---|
| 458 | } |
---|
| 459 | } |
---|
| 460 | else |
---|
| 461 | theGcd= nInit( 0 ); |
---|
| 462 | return theGcd; |
---|
| 463 | } |
---|
| 464 | |
---|
| 465 | number |
---|
| 466 | fglmVector::clearDenom() |
---|
| 467 | { |
---|
| 468 | number theLcm = nInit( 1 ); |
---|
| 469 | number current; |
---|
| 470 | BOOLEAN isZero = TRUE; |
---|
| 471 | int i; |
---|
| 472 | for ( i= size(); i > 0; i-- ) { |
---|
| 473 | if ( ! nIsZero( rep->getconstelem(i) ) ) { |
---|
| 474 | isZero= FALSE; |
---|
| 475 | number temp= nLcm( theLcm, rep->getconstelem( i ), currRing ); |
---|
| 476 | nDelete( &theLcm ); |
---|
| 477 | theLcm= temp; |
---|
| 478 | } |
---|
| 479 | } |
---|
| 480 | if ( isZero ) { |
---|
| 481 | nDelete( &theLcm ); |
---|
| 482 | theLcm= nInit( 0 ); |
---|
| 483 | } |
---|
| 484 | else { |
---|
| 485 | if ( ! nIsOne( theLcm ) ) { |
---|
| 486 | *this *= theLcm; |
---|
| 487 | for ( i= size(); i > 0; i-- ) { |
---|
| 488 | nNormalize( rep->getelem( i ) ); |
---|
| 489 | } |
---|
| 490 | } |
---|
| 491 | } |
---|
| 492 | return theLcm; |
---|
| 493 | } |
---|
| 494 | |
---|
| 495 | #endif |
---|
| 496 | // ---------------------------------------------------------------------------- |
---|
| 497 | // Local Variables: *** |
---|
| 498 | // compile-command: "make Singular" *** |
---|
| 499 | // page-delimiter: "^\\(\\|//!\\)" *** |
---|
| 500 | // fold-internal-margins: nil *** |
---|
| 501 | // End: *** |
---|