[7885020] | 1 | // ---------------------------------------------------------------------------- |
---|
| 2 | // spectrum.cc |
---|
| 3 | // begin of file |
---|
| 4 | // Stephan Endrass, endrass@mathematik.uni-mainz.de |
---|
| 5 | // 23.7.99 |
---|
| 6 | // ---------------------------------------------------------------------------- |
---|
| 7 | |
---|
| 8 | #define SPECTRUM_CC |
---|
| 9 | |
---|
| 10 | #include "mod2.h" |
---|
| 11 | |
---|
| 12 | #ifdef HAVE_SPECTRUM |
---|
| 13 | |
---|
| 14 | #ifdef SPECTRUM_PRINT |
---|
| 15 | #include <iostream.h> |
---|
| 16 | #ifndef SPECTRUM_IOSTREAM |
---|
| 17 | #include <stdio.h> |
---|
| 18 | #endif |
---|
| 19 | #endif |
---|
| 20 | |
---|
| 21 | #include <limits.h> |
---|
| 22 | |
---|
| 23 | #include "numbers.h" |
---|
| 24 | #include "polys.h" |
---|
| 25 | #include "ipid.h" |
---|
| 26 | #include "ideals.h" |
---|
| 27 | #include "kstd1.h" |
---|
| 28 | #include "stairc.h" |
---|
| 29 | #include "lists.h" |
---|
| 30 | #include "intvec.h" |
---|
| 31 | #include "ring.h" |
---|
| 32 | |
---|
| 33 | #include "multicnt.h" |
---|
| 34 | #include "GMPrat.h" |
---|
| 35 | #include "kmatrix.h" |
---|
| 36 | #include "npolygon.h" |
---|
| 37 | #include "splist.h" |
---|
| 38 | #include "semic.h" |
---|
| 39 | |
---|
| 40 | // ---------------------------------------------------------------------------- |
---|
[efaa52] | 41 | // test if the polynomial h has a term of total degree d |
---|
[7885020] | 42 | // ---------------------------------------------------------------------------- |
---|
| 43 | |
---|
[efaa52] | 44 | BOOLEAN hasTermOfDegree( poly h, int d ) |
---|
[7885020] | 45 | { |
---|
[efaa52] | 46 | do |
---|
| 47 | { |
---|
| 48 | if( pTotaldegree( h )== d ) |
---|
| 49 | return TRUE; |
---|
| 50 | pIter(h); |
---|
| 51 | } |
---|
| 52 | while( h!=(poly)NULL ); |
---|
| 53 | |
---|
| 54 | return FALSE; |
---|
| 55 | } |
---|
[7885020] | 56 | |
---|
[efaa52] | 57 | // ---------------------------------------------------------------------------- |
---|
| 58 | // test if the polynomial h has a constant term |
---|
| 59 | // ---------------------------------------------------------------------------- |
---|
[7885020] | 60 | |
---|
[efaa52] | 61 | static BOOLEAN inline hasConstTerm( poly h ) |
---|
| 62 | { |
---|
| 63 | return hasTermOfDegree(h,0); |
---|
[7885020] | 64 | } |
---|
| 65 | |
---|
| 66 | // ---------------------------------------------------------------------------- |
---|
| 67 | // test if the polynomial h has a linear term |
---|
| 68 | // ---------------------------------------------------------------------------- |
---|
| 69 | |
---|
[efaa52] | 70 | static BOOLEAN inline hasLinearTerm( poly h ) |
---|
[7885020] | 71 | { |
---|
[efaa52] | 72 | return hasTermOfDegree(h,1); |
---|
[7885020] | 73 | } |
---|
| 74 | |
---|
| 75 | // ---------------------------------------------------------------------------- |
---|
| 76 | // test if the ideal J has a lead monomial on the axis number k |
---|
| 77 | // ---------------------------------------------------------------------------- |
---|
| 78 | |
---|
| 79 | BOOLEAN hasAxis( ideal J,int k ) |
---|
| 80 | { |
---|
[efaa52] | 81 | int i; |
---|
[7885020] | 82 | |
---|
[efaa52] | 83 | for( i=0; i<IDELEMS(J); i++ ) |
---|
| 84 | { |
---|
| 85 | if (pIsPurePower( J->m[i] ) == k) return TRUE; |
---|
| 86 | } |
---|
| 87 | return FALSE; |
---|
[7885020] | 88 | } |
---|
| 89 | |
---|
| 90 | // ---------------------------------------------------------------------------- |
---|
| 91 | // test if the ideal J has 1 as a lead monomial |
---|
| 92 | // ---------------------------------------------------------------------------- |
---|
| 93 | |
---|
| 94 | int hasOne( ideal J ) |
---|
| 95 | { |
---|
[efaa52] | 96 | int i; |
---|
[7885020] | 97 | |
---|
[efaa52] | 98 | for( i=0; i<IDELEMS(J); i++ ) |
---|
| 99 | { |
---|
| 100 | if (pIsConstant(J->m[i])) return TRUE; |
---|
| 101 | } |
---|
| 102 | return FALSE; |
---|
[7885020] | 103 | } |
---|
| 104 | // ---------------------------------------------------------------------------- |
---|
| 105 | // test if m is a multiple of one of the monomials of f |
---|
| 106 | // ---------------------------------------------------------------------------- |
---|
| 107 | |
---|
| 108 | int isMultiple( poly f,poly m ) |
---|
| 109 | { |
---|
[efaa52] | 110 | while( f!=(poly)NULL ) |
---|
| 111 | { |
---|
| 112 | // --------------------------------------------------- |
---|
| 113 | // for a local order f|m is only possible if f>=m |
---|
| 114 | // --------------------------------------------------- |
---|
[7885020] | 115 | |
---|
[a6a239] | 116 | if( pLmCmp( f,m )>=0 ) |
---|
[efaa52] | 117 | { |
---|
[512a2b] | 118 | if( pLmDivisibleByNoComp( f,m ) ) |
---|
[efaa52] | 119 | { |
---|
| 120 | return TRUE; |
---|
| 121 | } |
---|
| 122 | else |
---|
| 123 | { |
---|
| 124 | pIter( f ); |
---|
| 125 | } |
---|
| 126 | } |
---|
| 127 | else |
---|
| 128 | { |
---|
| 129 | return FALSE; |
---|
[7885020] | 130 | } |
---|
[efaa52] | 131 | } |
---|
[7885020] | 132 | |
---|
[efaa52] | 133 | return FALSE; |
---|
[7885020] | 134 | } |
---|
| 135 | |
---|
| 136 | // ---------------------------------------------------------------------------- |
---|
| 137 | // compute the minimal monomial of minimmal weight>=max_weight |
---|
| 138 | // ---------------------------------------------------------------------------- |
---|
| 139 | |
---|
| 140 | poly computeWC( const newtonPolygon &np,Rational max_weight ) |
---|
| 141 | { |
---|
[efaa52] | 142 | poly m = pOne(); |
---|
| 143 | poly wc = (poly)NULL; |
---|
| 144 | int mdegree; |
---|
[7885020] | 145 | |
---|
[efaa52] | 146 | for( int i=1; i<=pVariables; i++ ) |
---|
| 147 | { |
---|
| 148 | mdegree = 1; |
---|
| 149 | pSetExp( m,i,mdegree ); |
---|
| 150 | // pSetm( m ); |
---|
| 151 | // np.weight_shift does not need pSetm( m ), postpone it |
---|
[7885020] | 152 | |
---|
[efaa52] | 153 | while( np.weight_shift( m )<max_weight ) |
---|
[7885020] | 154 | { |
---|
[efaa52] | 155 | mdegree++; |
---|
| 156 | pSetExp( m,i,mdegree ); |
---|
| 157 | // pSetm( m ); |
---|
| 158 | // np.weight_shift does not need pSetm( m ), postpone it |
---|
| 159 | } |
---|
| 160 | pSetm( m ); |
---|
[7885020] | 161 | |
---|
[a6a239] | 162 | if( i==1 || pCmp( m,wc )<0 ) |
---|
[efaa52] | 163 | { |
---|
| 164 | pDelete( &wc ); |
---|
| 165 | wc = pHead( m ); |
---|
[7885020] | 166 | } |
---|
| 167 | |
---|
[efaa52] | 168 | pSetExp( m,i,0 ); |
---|
| 169 | } |
---|
[7885020] | 170 | |
---|
[efaa52] | 171 | pDelete( &m ); |
---|
| 172 | |
---|
| 173 | return wc; |
---|
[7885020] | 174 | } |
---|
| 175 | |
---|
| 176 | // ---------------------------------------------------------------------------- |
---|
| 177 | // deletes all monomials of f which are less than hc |
---|
| 178 | // ---------------------------------------------------------------------------- |
---|
| 179 | |
---|
[efaa52] | 180 | static inline poly normalFormHC( poly f,poly hc ) |
---|
[7885020] | 181 | { |
---|
[efaa52] | 182 | poly *ptr = &f; |
---|
[7885020] | 183 | |
---|
[efaa52] | 184 | while( (*ptr)!=(poly)NULL ) |
---|
| 185 | { |
---|
[a6a239] | 186 | if( pLmCmp( *ptr,hc )>=0 ) |
---|
[7885020] | 187 | { |
---|
[efaa52] | 188 | ptr = &(pNext( *ptr )); |
---|
[7885020] | 189 | } |
---|
[efaa52] | 190 | else |
---|
| 191 | { |
---|
| 192 | pDelete( ptr ); |
---|
| 193 | return f; |
---|
| 194 | } |
---|
| 195 | } |
---|
[7885020] | 196 | |
---|
[efaa52] | 197 | return f; |
---|
[7885020] | 198 | } |
---|
| 199 | |
---|
| 200 | // ---------------------------------------------------------------------------- |
---|
| 201 | // deletes all monomials of f which are multiples of monomials of Z |
---|
| 202 | // ---------------------------------------------------------------------------- |
---|
| 203 | |
---|
[efaa52] | 204 | static inline poly normalFormZ( poly f,poly Z ) |
---|
[7885020] | 205 | { |
---|
[efaa52] | 206 | poly *ptr = &f; |
---|
[7885020] | 207 | |
---|
[efaa52] | 208 | while( (*ptr)!=(poly)NULL ) |
---|
| 209 | { |
---|
| 210 | if( !isMultiple( Z,*ptr ) ) |
---|
[7885020] | 211 | { |
---|
[efaa52] | 212 | ptr = &(pNext( *ptr )); |
---|
| 213 | } |
---|
| 214 | else |
---|
| 215 | { |
---|
[512a2b] | 216 | pDeleteLm(ptr); |
---|
[7885020] | 217 | } |
---|
[efaa52] | 218 | } |
---|
[7885020] | 219 | |
---|
[efaa52] | 220 | return f; |
---|
[7885020] | 221 | } |
---|
| 222 | |
---|
[efaa52] | 223 | // ?? HS: |
---|
| 224 | // Looks for the shortest polynomial f in stdJ which is divisible |
---|
| 225 | // by the monimial m, return its index in stdJ |
---|
[7885020] | 226 | // ---------------------------------------------------------------------------- |
---|
| 227 | // Looks for the first polynomial f in stdJ which satisfies m=LT(f)*eta |
---|
| 228 | // for a monomial eta. The return eta*f-m and cancel all monomials |
---|
| 229 | // which are smaller than the highest corner hc |
---|
| 230 | // ---------------------------------------------------------------------------- |
---|
| 231 | |
---|
[efaa52] | 232 | static inline int isLeadMonomial( poly m,ideal stdJ ) |
---|
[7885020] | 233 | { |
---|
[efaa52] | 234 | int length = INT_MAX; |
---|
| 235 | int result = -1; |
---|
[7885020] | 236 | |
---|
[efaa52] | 237 | for( int i=0; i<IDELEMS(stdJ); i++ ) |
---|
| 238 | { |
---|
[a6a239] | 239 | if( pCmp( stdJ->m[i],m )>=0 && pDivisibleBy( stdJ->m[i],m ) ) |
---|
[7885020] | 240 | { |
---|
[efaa52] | 241 | int tmp = pLength( stdJ->m[i] ); |
---|
[7885020] | 242 | |
---|
[efaa52] | 243 | if( tmp < length ) |
---|
| 244 | { |
---|
| 245 | length = tmp; |
---|
| 246 | result = i; |
---|
| 247 | } |
---|
[7885020] | 248 | } |
---|
[efaa52] | 249 | } |
---|
[7885020] | 250 | |
---|
[efaa52] | 251 | return result; |
---|
[7885020] | 252 | } |
---|
| 253 | |
---|
| 254 | // ---------------------------------------------------------------------------- |
---|
| 255 | // set the exponent of a monomial t an integer array |
---|
| 256 | // ---------------------------------------------------------------------------- |
---|
| 257 | |
---|
[efaa52] | 258 | static void setExp( poly m,int *r ) |
---|
[7885020] | 259 | { |
---|
[efaa52] | 260 | for( int i=pVariables; i>0; i-- ) |
---|
| 261 | { |
---|
| 262 | pSetExp( m,i,r[i-1] ); |
---|
| 263 | } |
---|
| 264 | pSetm( m ); |
---|
[7885020] | 265 | } |
---|
| 266 | |
---|
| 267 | // ---------------------------------------------------------------------------- |
---|
| 268 | // check if the ordering is a reverse wellordering, i.e. every monomial |
---|
| 269 | // is smaller than only finitely monomials |
---|
| 270 | // ---------------------------------------------------------------------------- |
---|
| 271 | |
---|
[efaa52] | 272 | static BOOLEAN isWell( void ) |
---|
[7885020] | 273 | { |
---|
[efaa52] | 274 | int b = rBlocks( currRing ); |
---|
| 275 | |
---|
| 276 | if( b==3 && |
---|
| 277 | ( currRing->order[0] == ringorder_ds || |
---|
| 278 | currRing->order[0] == ringorder_Ds || |
---|
| 279 | currRing->order[0] == ringorder_ws || |
---|
| 280 | currRing->order[0] == ringorder_Ws ) ) |
---|
| 281 | { |
---|
| 282 | return TRUE; |
---|
| 283 | } |
---|
[a7f305] | 284 | else if( b>=3 |
---|
| 285 | && (( currRing->order[0] ==ringorder_a |
---|
[efaa52] | 286 | && currRing->block1[0]==pVariables ) |
---|
[a7f305] | 287 | || (currRing->order[0]==ringorder_M |
---|
[efaa52] | 288 | && currRing->block1[0]==pVariables*pVariables ))) |
---|
| 289 | { |
---|
| 290 | for( int i=pVariables-1; i>=0; i-- ) |
---|
| 291 | { |
---|
| 292 | if( currRing->wvhdl[0][i]>=0 ) |
---|
| 293 | { |
---|
| 294 | return FALSE; |
---|
| 295 | } |
---|
[7885020] | 296 | } |
---|
[efaa52] | 297 | return TRUE; |
---|
| 298 | } |
---|
[7885020] | 299 | |
---|
[efaa52] | 300 | return FALSE; |
---|
[7885020] | 301 | } |
---|
| 302 | |
---|
| 303 | // ---------------------------------------------------------------------------- |
---|
| 304 | // compute all monomials not in stdJ and their normal forms |
---|
| 305 | // ---------------------------------------------------------------------------- |
---|
| 306 | |
---|
[efaa52] | 307 | static void computeNF( ideal stdJ,poly hc,poly wc,spectrumPolyList *NF ) |
---|
[7885020] | 308 | { |
---|
[efaa52] | 309 | int carry,k; |
---|
| 310 | multiCnt C( pVariables,0 ); |
---|
| 311 | poly Z = (poly)NULL; |
---|
[7885020] | 312 | |
---|
[efaa52] | 313 | int well = isWell( ); |
---|
[7885020] | 314 | |
---|
[efaa52] | 315 | do |
---|
| 316 | { |
---|
| 317 | poly m = pOne(); |
---|
| 318 | setExp( m,C.cnt ); |
---|
[7885020] | 319 | |
---|
[efaa52] | 320 | carry = FALSE; |
---|
[7885020] | 321 | |
---|
[efaa52] | 322 | k = isLeadMonomial( m,stdJ ); |
---|
[7885020] | 323 | |
---|
[efaa52] | 324 | if( k < 0 ) |
---|
| 325 | { |
---|
| 326 | // --------------------------- |
---|
| 327 | // m is not a lead monomial |
---|
| 328 | // --------------------------- |
---|
[7885020] | 329 | |
---|
[efaa52] | 330 | NF->insert_node( m,(poly)NULL ); |
---|
| 331 | } |
---|
| 332 | else if( isMultiple( Z,m ) ) |
---|
| 333 | { |
---|
| 334 | // ------------------------------------ |
---|
| 335 | // m is trivially in the ideal stdJ |
---|
| 336 | // ------------------------------------ |
---|
[7885020] | 337 | |
---|
[efaa52] | 338 | pDelete( &m ); |
---|
| 339 | carry = TRUE; |
---|
| 340 | } |
---|
[a6a239] | 341 | else if( pCmp( m,hc ) < 0 || pCmp( m,wc ) < 0 ) |
---|
[efaa52] | 342 | { |
---|
| 343 | // ------------------- |
---|
| 344 | // we do not need m |
---|
| 345 | // ------------------- |
---|
[7885020] | 346 | |
---|
[efaa52] | 347 | pDelete( &m ); |
---|
| 348 | carry = TRUE; |
---|
| 349 | } |
---|
| 350 | else |
---|
| 351 | { |
---|
| 352 | // -------------------------- |
---|
| 353 | // compute lazy normal form |
---|
| 354 | // -------------------------- |
---|
| 355 | |
---|
| 356 | poly multiplicant = pDivide( m,stdJ->m[k] ); |
---|
| 357 | pGetCoeff( multiplicant ) = nInit(1); |
---|
[7885020] | 358 | |
---|
[a6a239] | 359 | poly nf = pMult_mm( pCopy( stdJ->m[k] ), multiplicant ); |
---|
[7885020] | 360 | |
---|
[efaa52] | 361 | pDelete( &multiplicant ); |
---|
[7885020] | 362 | |
---|
[efaa52] | 363 | nf = normalFormHC( nf,hc ); |
---|
[7885020] | 364 | |
---|
[efaa52] | 365 | if( pNext( nf )==(poly)NULL ) |
---|
| 366 | { |
---|
| 367 | // ---------------------------------- |
---|
| 368 | // normal form of m is m itself |
---|
| 369 | // ---------------------------------- |
---|
[7885020] | 370 | |
---|
[efaa52] | 371 | pDelete( &nf ); |
---|
| 372 | NF->delete_monomial( m ); |
---|
| 373 | Z = pAdd( Z,m ); |
---|
| 374 | carry = TRUE; |
---|
| 375 | } |
---|
| 376 | else |
---|
| 377 | { |
---|
| 378 | nf = normalFormZ( nf,Z ); |
---|
| 379 | |
---|
| 380 | if( pNext( nf )==(poly)NULL ) |
---|
| 381 | { |
---|
| 382 | // ---------------------------------- |
---|
| 383 | // normal form of m is m itself |
---|
| 384 | // ---------------------------------- |
---|
| 385 | |
---|
| 386 | pDelete( &nf ); |
---|
| 387 | NF->delete_monomial( m ); |
---|
| 388 | Z = pAdd( Z,m ); |
---|
| 389 | carry = TRUE; |
---|
| 390 | } |
---|
| 391 | else |
---|
| 392 | { |
---|
| 393 | // ------------------------------------ |
---|
| 394 | // normal form of m is a polynomial |
---|
| 395 | // ------------------------------------ |
---|
| 396 | |
---|
| 397 | pNorm( nf ); |
---|
| 398 | if( well==TRUE ) |
---|
| 399 | { |
---|
| 400 | NF->insert_node( m,nf ); |
---|
| 401 | } |
---|
| 402 | else |
---|
| 403 | { |
---|
| 404 | poly nfhard = kNF( stdJ,(ideal)NULL,pNext( nf ),0,0 ); |
---|
| 405 | nfhard = normalFormHC( nfhard,hc ); |
---|
| 406 | nfhard = normalFormZ ( nfhard,Z ); |
---|
| 407 | |
---|
| 408 | if( nfhard==(poly)NULL ) |
---|
[7885020] | 409 | { |
---|
[efaa52] | 410 | NF->delete_monomial( m ); |
---|
| 411 | Z = pAdd( Z,m ); |
---|
| 412 | carry = TRUE; |
---|
[7885020] | 413 | } |
---|
| 414 | else |
---|
| 415 | { |
---|
[efaa52] | 416 | pDelete( &pNext( nf ) ); |
---|
| 417 | pNext( nf ) = nfhard; |
---|
| 418 | NF->insert_node( m,nf ); |
---|
[7885020] | 419 | } |
---|
[efaa52] | 420 | } |
---|
[7885020] | 421 | } |
---|
[efaa52] | 422 | } |
---|
[7885020] | 423 | } |
---|
[efaa52] | 424 | } |
---|
| 425 | while( C.inc( carry ) ); |
---|
[7885020] | 426 | |
---|
[efaa52] | 427 | // delete single monomials |
---|
[7885020] | 428 | |
---|
[efaa52] | 429 | BOOLEAN not_finished; |
---|
[7885020] | 430 | |
---|
[efaa52] | 431 | do |
---|
| 432 | { |
---|
| 433 | not_finished = FALSE; |
---|
[7885020] | 434 | |
---|
[efaa52] | 435 | spectrumPolyNode *node = NF->root; |
---|
[7885020] | 436 | |
---|
[efaa52] | 437 | while( node!=(spectrumPolyNode*)NULL ) |
---|
| 438 | { |
---|
| 439 | if( node->nf!=(poly)NULL && pNext( node->nf )==(poly)NULL ) |
---|
| 440 | { |
---|
| 441 | NF->delete_monomial( node->nf ); |
---|
| 442 | not_finished = TRUE; |
---|
| 443 | node = (spectrumPolyNode*)NULL; |
---|
| 444 | } |
---|
| 445 | else |
---|
| 446 | { |
---|
| 447 | node = node->next; |
---|
| 448 | } |
---|
[7885020] | 449 | } |
---|
[efaa52] | 450 | } while( not_finished ); |
---|
[7885020] | 451 | |
---|
[efaa52] | 452 | pDelete( &Z ); |
---|
[7885020] | 453 | } |
---|
| 454 | |
---|
| 455 | // ---------------------------------------------------------------------------- |
---|
| 456 | // this is the main spectrum computation function |
---|
| 457 | // ---------------------------------------------------------------------------- |
---|
| 458 | |
---|
[efaa52] | 459 | static spectrumState spectrumCompute( poly h,lists *L,int fast ) |
---|
[7885020] | 460 | { |
---|
[efaa52] | 461 | int i,j; |
---|
| 462 | |
---|
| 463 | #ifdef SPECTRUM_DEBUG |
---|
| 464 | #ifdef SPECTRUM_PRINT |
---|
| 465 | #ifdef SPECTRUM_IOSTREAM |
---|
| 466 | cout << "spectrumCompute\n"; |
---|
| 467 | if( fast==0 ) cout << " no optimization" << endl; |
---|
| 468 | if( fast==1 ) cout << " weight optimization" << endl; |
---|
| 469 | if( fast==2 ) cout << " symmetry optimization" << endl; |
---|
| 470 | #else |
---|
| 471 | fprintf( stdout,"spectrumCompute\n" ); |
---|
| 472 | if( fast==0 ) fprintf( stdout," no optimization\n" ); |
---|
| 473 | if( fast==1 ) fprintf( stdout," weight optimization\n" ); |
---|
| 474 | if( fast==2 ) fprintf( stdout," symmetry optimization\n" ); |
---|
| 475 | #endif |
---|
| 476 | #endif |
---|
| 477 | #endif |
---|
| 478 | |
---|
| 479 | // ---------------------- |
---|
| 480 | // check if h is zero |
---|
| 481 | // ---------------------- |
---|
| 482 | |
---|
| 483 | if( h==(poly)NULL ) |
---|
| 484 | { |
---|
| 485 | return spectrumZero; |
---|
| 486 | } |
---|
| 487 | |
---|
| 488 | // ---------------------------------- |
---|
| 489 | // check if h has a constant term |
---|
| 490 | // ---------------------------------- |
---|
| 491 | |
---|
| 492 | if( hasConstTerm( h ) ) |
---|
| 493 | { |
---|
| 494 | return spectrumBadPoly; |
---|
| 495 | } |
---|
| 496 | |
---|
| 497 | // -------------------------------- |
---|
| 498 | // check if h has a linear term |
---|
| 499 | // -------------------------------- |
---|
| 500 | |
---|
| 501 | if( hasLinearTerm( h ) ) |
---|
| 502 | { |
---|
[c232af] | 503 | *L = (lists)omAllocBin( slists _bin); |
---|
[efaa52] | 504 | (*L)->Init( 1 ); |
---|
| 505 | (*L)->m[0].rtyp = INT_CMD; // milnor number |
---|
| 506 | /* (*L)->m[0].data = (void*)0;a -- done by Init */ |
---|
| 507 | |
---|
| 508 | return spectrumNoSingularity; |
---|
| 509 | } |
---|
| 510 | |
---|
| 511 | // ---------------------------------- |
---|
| 512 | // compute the jacobi ideal of (h) |
---|
| 513 | // ---------------------------------- |
---|
| 514 | |
---|
| 515 | ideal J = NULL; |
---|
| 516 | J = idInit( pVariables,1 ); |
---|
| 517 | |
---|
| 518 | #ifdef SPECTRUM_DEBUG |
---|
| 519 | #ifdef SPECTRUM_PRINT |
---|
| 520 | #ifdef SPECTRUM_IOSTREAM |
---|
| 521 | cout << "\n computing the Jacobi ideal...\n"; |
---|
| 522 | #else |
---|
| 523 | fprintf( stdout,"\n computing the Jacobi ideal...\n" ); |
---|
| 524 | #endif |
---|
| 525 | #endif |
---|
| 526 | #endif |
---|
| 527 | |
---|
| 528 | for( i=0; i<pVariables; i++ ) |
---|
| 529 | { |
---|
| 530 | J->m[i] = pDiff( h,i+1); //j ); |
---|
[7885020] | 531 | |
---|
| 532 | #ifdef SPECTRUM_DEBUG |
---|
| 533 | #ifdef SPECTRUM_PRINT |
---|
| 534 | #ifdef SPECTRUM_IOSTREAM |
---|
[efaa52] | 535 | cout << " "; |
---|
[7885020] | 536 | #else |
---|
[efaa52] | 537 | fprintf( stdout," " ); |
---|
[7885020] | 538 | #endif |
---|
[efaa52] | 539 | pWrite( J->m[i] ); |
---|
[7885020] | 540 | #endif |
---|
| 541 | #endif |
---|
[efaa52] | 542 | } |
---|
| 543 | |
---|
| 544 | // -------------------------------------------- |
---|
| 545 | // compute a standard basis stdJ of jac(h) |
---|
| 546 | // -------------------------------------------- |
---|
| 547 | |
---|
| 548 | #ifdef SPECTRUM_DEBUG |
---|
| 549 | #ifdef SPECTRUM_PRINT |
---|
| 550 | #ifdef SPECTRUM_IOSTREAM |
---|
| 551 | cout << endl; |
---|
| 552 | cout << " computing a standard basis..." << endl; |
---|
| 553 | #else |
---|
| 554 | fprintf( stdout,"\n" ); |
---|
| 555 | fprintf( stdout," computing a standard basis...\n" ); |
---|
| 556 | #endif |
---|
| 557 | #endif |
---|
| 558 | #endif |
---|
| 559 | |
---|
| 560 | ideal stdJ = kStd(J,currQuotient,isNotHomog,NULL); |
---|
| 561 | idSkipZeroes( stdJ ); |
---|
| 562 | |
---|
| 563 | #ifdef SPECTRUM_DEBUG |
---|
| 564 | #ifdef SPECTRUM_PRINT |
---|
| 565 | for( i=0; i<IDELEMS(stdJ); i++ ) |
---|
| 566 | { |
---|
| 567 | #ifdef SPECTRUM_IOSTREAM |
---|
| 568 | cout << " "; |
---|
| 569 | #else |
---|
| 570 | fprintf( stdout," " ); |
---|
| 571 | #endif |
---|
[7885020] | 572 | |
---|
[efaa52] | 573 | pWrite( stdJ->m[i] ); |
---|
[7885020] | 574 | } |
---|
[efaa52] | 575 | #endif |
---|
| 576 | #endif |
---|
[7885020] | 577 | |
---|
[efaa52] | 578 | idDelete( &J ); |
---|
[7885020] | 579 | |
---|
[efaa52] | 580 | // ------------------------------------------ |
---|
| 581 | // check if the h has a singularity |
---|
| 582 | // ------------------------------------------ |
---|
[7885020] | 583 | |
---|
[efaa52] | 584 | if( hasOne( stdJ ) ) |
---|
| 585 | { |
---|
| 586 | // ------------------------------- |
---|
| 587 | // h is smooth in the origin |
---|
| 588 | // return only the Milnor number |
---|
| 589 | // ------------------------------- |
---|
[7885020] | 590 | |
---|
[c232af] | 591 | *L = (lists)omAllocBin( slists _bin); |
---|
[efaa52] | 592 | (*L)->Init( 1 ); |
---|
| 593 | (*L)->m[0].rtyp = INT_CMD; // milnor number |
---|
| 594 | /* (*L)->m[0].data = (void*)0;a -- done by Init */ |
---|
[7885020] | 595 | |
---|
[efaa52] | 596 | return spectrumNoSingularity; |
---|
| 597 | } |
---|
[7885020] | 598 | |
---|
[efaa52] | 599 | // ------------------------------------------ |
---|
| 600 | // check if the singularity h is isolated |
---|
| 601 | // ------------------------------------------ |
---|
[7885020] | 602 | |
---|
[efaa52] | 603 | for( i=pVariables; i>0; i-- ) |
---|
| 604 | { |
---|
| 605 | if( hasAxis( stdJ,i )==FALSE ) |
---|
[7885020] | 606 | { |
---|
[efaa52] | 607 | return spectrumNotIsolated; |
---|
[7885020] | 608 | } |
---|
[efaa52] | 609 | } |
---|
[7885020] | 610 | |
---|
[efaa52] | 611 | // ------------------------------------------ |
---|
| 612 | // compute the highest corner hc of stdJ |
---|
| 613 | // ------------------------------------------ |
---|
[7885020] | 614 | |
---|
[efaa52] | 615 | #ifdef SPECTRUM_DEBUG |
---|
| 616 | #ifdef SPECTRUM_PRINT |
---|
| 617 | #ifdef SPECTRUM_IOSTREAM |
---|
| 618 | cout << "\n computing the highest corner...\n"; |
---|
| 619 | #else |
---|
| 620 | fprintf( stdout,"\n computing the highest corner...\n" ); |
---|
| 621 | #endif |
---|
| 622 | #endif |
---|
| 623 | #endif |
---|
[7885020] | 624 | |
---|
[efaa52] | 625 | poly hc = (poly)NULL; |
---|
[7885020] | 626 | |
---|
[efaa52] | 627 | scComputeHC( stdJ,0,hc ); |
---|
[7885020] | 628 | |
---|
[efaa52] | 629 | if( hc!=(poly)NULL ) |
---|
| 630 | { |
---|
| 631 | pGetCoeff(hc) = nInit(1); |
---|
[7885020] | 632 | |
---|
| 633 | for( i=pVariables; i>0; i-- ) |
---|
| 634 | { |
---|
[efaa52] | 635 | if( pGetExp( hc,i )>0 ) pDecrExp( hc,i ); |
---|
| 636 | } |
---|
| 637 | pSetm( hc ); |
---|
| 638 | } |
---|
| 639 | else |
---|
| 640 | { |
---|
| 641 | return spectrumNoHC; |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | #ifdef SPECTRUM_DEBUG |
---|
| 645 | #ifdef SPECTRUM_PRINT |
---|
| 646 | #ifdef SPECTRUM_IOSTREAM |
---|
| 647 | cout << " "; |
---|
| 648 | #else |
---|
| 649 | fprintf( stdout," " ); |
---|
| 650 | #endif |
---|
| 651 | pWrite( hc ); |
---|
| 652 | #endif |
---|
| 653 | #endif |
---|
| 654 | |
---|
| 655 | // ---------------------------------------- |
---|
| 656 | // compute the Newton polygon nph of h |
---|
| 657 | // ---------------------------------------- |
---|
| 658 | |
---|
| 659 | #ifdef SPECTRUM_DEBUG |
---|
| 660 | #ifdef SPECTRUM_PRINT |
---|
| 661 | #ifdef SPECTRUM_IOSTREAM |
---|
| 662 | cout << "\n computing the newton polygon...\n"; |
---|
| 663 | #else |
---|
| 664 | fprintf( stdout,"\n computing the newton polygon...\n" ); |
---|
| 665 | #endif |
---|
| 666 | #endif |
---|
| 667 | #endif |
---|
| 668 | |
---|
| 669 | newtonPolygon nph( h ); |
---|
| 670 | |
---|
| 671 | #ifdef SPECTRUM_DEBUG |
---|
| 672 | #ifdef SPECTRUM_PRINT |
---|
| 673 | cout << nph; |
---|
| 674 | #endif |
---|
| 675 | #endif |
---|
| 676 | |
---|
| 677 | // ----------------------------------------------- |
---|
| 678 | // compute the weight corner wc of (stdj,nph) |
---|
| 679 | // ----------------------------------------------- |
---|
| 680 | |
---|
| 681 | #ifdef SPECTRUM_DEBUG |
---|
| 682 | #ifdef SPECTRUM_PRINT |
---|
| 683 | #ifdef SPECTRUM_IOSTREAM |
---|
| 684 | cout << "\n computing the weight corner...\n"; |
---|
| 685 | #else |
---|
| 686 | fprintf( stdout,"\n computing the weight corner...\n" ); |
---|
| 687 | #endif |
---|
| 688 | #endif |
---|
| 689 | #endif |
---|
| 690 | |
---|
| 691 | poly wc = ( fast==0 ? pCopy( hc ) : |
---|
| 692 | ( fast==1 ? computeWC( nph,(Rational)pVariables ) : |
---|
| 693 | /* fast==2 */computeWC( nph,((Rational)pVariables)/(Rational)2 ) ) ); |
---|
| 694 | |
---|
| 695 | #ifdef SPECTRUM_DEBUG |
---|
| 696 | #ifdef SPECTRUM_PRINT |
---|
| 697 | #ifdef SPECTRUM_IOSTREAM |
---|
| 698 | cout << " "; |
---|
| 699 | #else |
---|
| 700 | fprintf( stdout," " ); |
---|
| 701 | #endif |
---|
| 702 | pWrite( wc ); |
---|
| 703 | #endif |
---|
| 704 | #endif |
---|
| 705 | |
---|
| 706 | // ------------- |
---|
| 707 | // compute NF |
---|
| 708 | // ------------- |
---|
| 709 | |
---|
| 710 | #ifdef SPECTRUM_DEBUG |
---|
| 711 | #ifdef SPECTRUM_PRINT |
---|
| 712 | #ifdef SPECTRUM_IOSTREAM |
---|
| 713 | cout << "\n computing NF...\n" << endl; |
---|
| 714 | #else |
---|
| 715 | fprintf( stdout,"\n computing NF...\n" ); |
---|
| 716 | #endif |
---|
| 717 | #endif |
---|
| 718 | #endif |
---|
| 719 | |
---|
| 720 | spectrumPolyList NF( &nph ); |
---|
| 721 | |
---|
| 722 | computeNF( stdJ,hc,wc,&NF ); |
---|
| 723 | |
---|
| 724 | #ifdef SPECTRUM_DEBUG |
---|
| 725 | #ifdef SPECTRUM_PRINT |
---|
| 726 | cout << NF; |
---|
| 727 | #ifdef SPECTRUM_IOSTREAM |
---|
| 728 | cout << endl; |
---|
| 729 | #else |
---|
| 730 | fprintf( stdout,"\n" ); |
---|
| 731 | #endif |
---|
| 732 | #endif |
---|
| 733 | #endif |
---|
| 734 | |
---|
| 735 | // ---------------------------- |
---|
| 736 | // compute the spectrum of h |
---|
| 737 | // ---------------------------- |
---|
| 738 | |
---|
| 739 | return NF.spectrum( L,fast ); |
---|
[7885020] | 740 | } |
---|
| 741 | |
---|
| 742 | // ---------------------------------------------------------------------------- |
---|
| 743 | // check if currRing is local |
---|
| 744 | // ---------------------------------------------------------------------------- |
---|
| 745 | |
---|
[efaa52] | 746 | static BOOLEAN ringIsLocal( void ) |
---|
[7885020] | 747 | { |
---|
[efaa52] | 748 | poly m = pOne(); |
---|
| 749 | poly one = pOne(); |
---|
[7885020] | 750 | |
---|
[efaa52] | 751 | for( int i=pVariables; i>0; i-- ) |
---|
| 752 | { |
---|
| 753 | pSetExp( m,i,1 ); |
---|
| 754 | pSetm( m ); |
---|
[7885020] | 755 | |
---|
[a6a239] | 756 | if( pCmp( m,one )>0 ) |
---|
[efaa52] | 757 | { |
---|
| 758 | return FALSE; |
---|
[7885020] | 759 | } |
---|
[efaa52] | 760 | pSetExp( m,i,0 ); |
---|
| 761 | } |
---|
[7885020] | 762 | |
---|
[efaa52] | 763 | pDelete( &m ); |
---|
| 764 | pDelete( &one ); |
---|
[7885020] | 765 | |
---|
[efaa52] | 766 | return TRUE; |
---|
[7885020] | 767 | } |
---|
| 768 | |
---|
[efaa52] | 769 | // ---------------------------------------------------------------------------- |
---|
| 770 | // print error message corresponding to spectrumState state: |
---|
| 771 | // ---------------------------------------------------------------------------- |
---|
| 772 | static void spectrumPrintError(spectrumState state) |
---|
| 773 | { |
---|
| 774 | switch( state ) |
---|
| 775 | { |
---|
| 776 | case spectrumZero: |
---|
| 777 | WerrorS( "polynomial is zero" ); |
---|
| 778 | break; |
---|
| 779 | case spectrumBadPoly: |
---|
| 780 | WerrorS( "polynomial has constant term" ); |
---|
| 781 | break; |
---|
| 782 | case spectrumNoSingularity: |
---|
| 783 | WerrorS( "not a singularity" ); |
---|
| 784 | break; |
---|
| 785 | case spectrumNotIsolated: |
---|
| 786 | WerrorS( "the singularity is not isolated" ); |
---|
| 787 | break; |
---|
| 788 | case spectrumNoHC: |
---|
| 789 | WerrorS( "highest corner cannot be computed" ); |
---|
| 790 | break; |
---|
| 791 | case spectrumDegenerate: |
---|
| 792 | WerrorS( "principal part is degenerate" ); |
---|
| 793 | break; |
---|
| 794 | case spectrumOK: |
---|
| 795 | break; |
---|
| 796 | |
---|
| 797 | default: |
---|
| 798 | WerrorS( "unknown error occurred" ); |
---|
| 799 | break; |
---|
| 800 | } |
---|
| 801 | } |
---|
[7885020] | 802 | // ---------------------------------------------------------------------------- |
---|
| 803 | // this procedure is called from the interpreter |
---|
| 804 | // ---------------------------------------------------------------------------- |
---|
| 805 | // first = polynomial |
---|
| 806 | // result = list of spectrum numbers |
---|
| 807 | // ---------------------------------------------------------------------------- |
---|
| 808 | |
---|
| 809 | BOOLEAN spectrumProc( leftv result,leftv first ) |
---|
| 810 | { |
---|
[efaa52] | 811 | spectrumState state = spectrumOK; |
---|
[7885020] | 812 | |
---|
[efaa52] | 813 | // ------------------- |
---|
| 814 | // check consistency |
---|
| 815 | // ------------------- |
---|
[7885020] | 816 | |
---|
[efaa52] | 817 | // check for a local ring |
---|
[7885020] | 818 | |
---|
[efaa52] | 819 | if( !ringIsLocal( ) ) |
---|
| 820 | { |
---|
| 821 | WerrorS( "only works for local orderings" ); |
---|
| 822 | state = spectrumWrongRing; |
---|
| 823 | } |
---|
| 824 | |
---|
| 825 | // no quotient rings are allowed |
---|
| 826 | |
---|
| 827 | else if( currRing->qideal != NULL ) |
---|
| 828 | { |
---|
| 829 | WerrorS( "does not work in quotient rings" ); |
---|
| 830 | state = spectrumWrongRing; |
---|
| 831 | } |
---|
| 832 | else |
---|
| 833 | { |
---|
| 834 | lists L = (lists)NULL; |
---|
| 835 | int flag = 1; // weight corner optimization is safe |
---|
[7885020] | 836 | |
---|
[efaa52] | 837 | state = spectrumCompute( (poly)first->Data( ),&L,flag ); |
---|
[7885020] | 838 | |
---|
[efaa52] | 839 | if( state==spectrumOK ) |
---|
[7885020] | 840 | { |
---|
[efaa52] | 841 | result->rtyp = LIST_CMD; |
---|
| 842 | result->data = (char*)L; |
---|
[7885020] | 843 | } |
---|
| 844 | else |
---|
| 845 | { |
---|
[efaa52] | 846 | spectrumPrintError(state); |
---|
[7885020] | 847 | } |
---|
[efaa52] | 848 | } |
---|
[7885020] | 849 | |
---|
[efaa52] | 850 | return (state!=spectrumOK); |
---|
[7885020] | 851 | } |
---|
| 852 | |
---|
| 853 | // ---------------------------------------------------------------------------- |
---|
| 854 | // this procedure is called from the interpreter |
---|
| 855 | // ---------------------------------------------------------------------------- |
---|
| 856 | // first = polynomial |
---|
| 857 | // result = list of spectrum numbers |
---|
| 858 | // ---------------------------------------------------------------------------- |
---|
| 859 | |
---|
| 860 | BOOLEAN spectrumfProc( leftv result,leftv first ) |
---|
| 861 | { |
---|
[efaa52] | 862 | spectrumState state = spectrumOK; |
---|
| 863 | |
---|
| 864 | // ------------------- |
---|
| 865 | // check consistency |
---|
| 866 | // ------------------- |
---|
| 867 | |
---|
| 868 | // check for a local polynomial ring |
---|
| 869 | |
---|
[a7f305] | 870 | if( currRing->OrdSgn != -1 ) |
---|
[efaa52] | 871 | // ?? HS: the test above is also true for k[x][[y]], k[[x]][y] |
---|
[5c2863] | 872 | // or should we use: |
---|
[efaa52] | 873 | //if( !ringIsLocal( ) ) |
---|
| 874 | { |
---|
| 875 | WerrorS( "only works for local orderings" ); |
---|
| 876 | state = spectrumWrongRing; |
---|
| 877 | } |
---|
| 878 | else if( currRing->qideal != NULL ) |
---|
| 879 | { |
---|
| 880 | WerrorS( "does not work in quotient rings" ); |
---|
| 881 | state = spectrumWrongRing; |
---|
| 882 | } |
---|
| 883 | else |
---|
| 884 | { |
---|
| 885 | lists L = (lists)NULL; |
---|
| 886 | int flag = 2; // symmetric optimization |
---|
| 887 | |
---|
| 888 | state = spectrumCompute( (poly)first->Data( ),&L,flag ); |
---|
| 889 | |
---|
| 890 | if( state==spectrumOK ) |
---|
| 891 | { |
---|
| 892 | result->rtyp = LIST_CMD; |
---|
| 893 | result->data = (char*)L; |
---|
[7885020] | 894 | } |
---|
| 895 | else |
---|
| 896 | { |
---|
[efaa52] | 897 | spectrumPrintError(state); |
---|
[7885020] | 898 | } |
---|
[efaa52] | 899 | } |
---|
[7885020] | 900 | |
---|
[efaa52] | 901 | return (state!=spectrumOK); |
---|
[7885020] | 902 | } |
---|
| 903 | |
---|
| 904 | // ---------------------------------------------------------------------------- |
---|
| 905 | // check if a list is a spectrum |
---|
| 906 | // check for: |
---|
| 907 | // list has 6 elements |
---|
| 908 | // 1st element is int (mu=Milnor number) |
---|
| 909 | // 2nd element is int (pg=geometrical genus) |
---|
| 910 | // 3rd element is int (n =number of different spectrum numbers) |
---|
| 911 | // 4th element is intvec (num=numerators) |
---|
| 912 | // 5th element is intvec (den=denomiantors) |
---|
| 913 | // 6th element is intvec (mul=multiplicities) |
---|
| 914 | // exactly n numerators |
---|
| 915 | // exactly n denominators |
---|
| 916 | // exactly n multiplicities |
---|
| 917 | // mu>0 |
---|
| 918 | // pg>=0 |
---|
| 919 | // n>0 |
---|
| 920 | // num>0 |
---|
| 921 | // den>0 |
---|
| 922 | // mul>0 |
---|
| 923 | // symmetriy with respect to numberofvariables/2 |
---|
| 924 | // monotony |
---|
| 925 | // mu = sum of all multiplicities |
---|
| 926 | // pg = sum of all multiplicities where num/den<=1 |
---|
| 927 | // ---------------------------------------------------------------------------- |
---|
| 928 | |
---|
| 929 | semicState list_is_spectrum( lists l ) |
---|
| 930 | { |
---|
| 931 | // ------------------- |
---|
| 932 | // check list length |
---|
| 933 | // ------------------- |
---|
| 934 | |
---|
| 935 | if( l->nr < 5 ) |
---|
| 936 | { |
---|
| 937 | return semicListTooShort; |
---|
| 938 | } |
---|
| 939 | else if( l->nr > 5 ) |
---|
| 940 | { |
---|
| 941 | return semicListTooLong; |
---|
| 942 | } |
---|
| 943 | |
---|
| 944 | // ------------- |
---|
| 945 | // check types |
---|
| 946 | // ------------- |
---|
| 947 | |
---|
| 948 | if( l->m[0].rtyp != INT_CMD ) |
---|
| 949 | { |
---|
| 950 | return semicListFirstElementWrongType; |
---|
| 951 | } |
---|
| 952 | else if( l->m[1].rtyp != INT_CMD ) |
---|
| 953 | { |
---|
| 954 | return semicListSecondElementWrongType; |
---|
| 955 | } |
---|
| 956 | else if( l->m[2].rtyp != INT_CMD ) |
---|
| 957 | { |
---|
| 958 | return semicListThirdElementWrongType; |
---|
| 959 | } |
---|
| 960 | else if( l->m[3].rtyp != INTVEC_CMD ) |
---|
| 961 | { |
---|
| 962 | return semicListFourthElementWrongType; |
---|
| 963 | } |
---|
| 964 | else if( l->m[4].rtyp != INTVEC_CMD ) |
---|
| 965 | { |
---|
| 966 | return semicListFifthElementWrongType; |
---|
| 967 | } |
---|
| 968 | else if( l->m[5].rtyp != INTVEC_CMD ) |
---|
| 969 | { |
---|
| 970 | return semicListSixthElementWrongType; |
---|
| 971 | } |
---|
| 972 | |
---|
| 973 | // ------------------------- |
---|
| 974 | // check number of entries |
---|
| 975 | // ------------------------- |
---|
| 976 | |
---|
| 977 | int mu = (int)(l->m[0].Data( )); |
---|
| 978 | int pg = (int)(l->m[1].Data( )); |
---|
| 979 | int n = (int)(l->m[2].Data( )); |
---|
| 980 | |
---|
| 981 | if( n <= 0 ) |
---|
| 982 | { |
---|
| 983 | return semicListNNegative; |
---|
| 984 | } |
---|
| 985 | |
---|
| 986 | intvec *num = (intvec*)l->m[3].Data( ); |
---|
| 987 | intvec *den = (intvec*)l->m[4].Data( ); |
---|
| 988 | intvec *mul = (intvec*)l->m[5].Data( ); |
---|
| 989 | |
---|
| 990 | if( n != num->length( ) ) |
---|
| 991 | { |
---|
| 992 | return semicListWrongNumberOfNumerators; |
---|
| 993 | } |
---|
| 994 | else if( n != den->length( ) ) |
---|
| 995 | { |
---|
| 996 | return semicListWrongNumberOfDenominators; |
---|
| 997 | } |
---|
| 998 | else if( n != mul->length( ) ) |
---|
| 999 | { |
---|
| 1000 | return semicListWrongNumberOfMultiplicities; |
---|
| 1001 | } |
---|
| 1002 | |
---|
| 1003 | // -------- |
---|
| 1004 | // values |
---|
| 1005 | // -------- |
---|
| 1006 | |
---|
| 1007 | if( mu <= 0 ) |
---|
| 1008 | { |
---|
| 1009 | return semicListMuNegative; |
---|
| 1010 | } |
---|
| 1011 | if( pg < 0 ) |
---|
| 1012 | { |
---|
| 1013 | return semicListPgNegative; |
---|
| 1014 | } |
---|
| 1015 | |
---|
| 1016 | int i; |
---|
| 1017 | |
---|
| 1018 | for( i=0; i<n; i++ ) |
---|
| 1019 | { |
---|
| 1020 | if( (*num)[i] <= 0 ) |
---|
| 1021 | { |
---|
| 1022 | return semicListNumNegative; |
---|
| 1023 | } |
---|
| 1024 | if( (*den)[i] <= 0 ) |
---|
| 1025 | { |
---|
| 1026 | return semicListDenNegative; |
---|
| 1027 | } |
---|
| 1028 | if( (*mul)[i] <= 0 ) |
---|
| 1029 | { |
---|
| 1030 | return semicListMulNegative; |
---|
| 1031 | } |
---|
| 1032 | } |
---|
| 1033 | |
---|
| 1034 | // ---------------- |
---|
| 1035 | // check symmetry |
---|
| 1036 | // ---------------- |
---|
| 1037 | |
---|
| 1038 | int j; |
---|
| 1039 | |
---|
| 1040 | for( i=0, j=n-1; i<=j; i++,j-- ) |
---|
| 1041 | { |
---|
| 1042 | if( (*num)[i] != pVariables*((*den)[i]) - (*num)[j] || |
---|
| 1043 | (*den)[i] != (*den)[j] || |
---|
| 1044 | (*mul)[i] != (*mul)[j] ) |
---|
| 1045 | { |
---|
| 1046 | return semicListNotSymmetric; |
---|
| 1047 | } |
---|
| 1048 | } |
---|
| 1049 | |
---|
| 1050 | // ---------------- |
---|
| 1051 | // check monotony |
---|
| 1052 | // ---------------- |
---|
| 1053 | |
---|
| 1054 | for( i=0, j=1; i<n/2; i++,j++ ) |
---|
| 1055 | { |
---|
| 1056 | if( (*num)[i]*(*den)[j] >= (*num)[j]*(*den)[i] ) |
---|
| 1057 | { |
---|
| 1058 | return semicListNotMonotonous; |
---|
| 1059 | } |
---|
| 1060 | } |
---|
| 1061 | |
---|
| 1062 | // --------------------- |
---|
| 1063 | // check Milnor number |
---|
| 1064 | // --------------------- |
---|
| 1065 | |
---|
| 1066 | for( mu=0, i=0; i<n; i++ ) |
---|
| 1067 | { |
---|
| 1068 | mu += (*mul)[i]; |
---|
| 1069 | } |
---|
| 1070 | |
---|
| 1071 | if( mu != (int)(l->m[0].Data( )) ) |
---|
| 1072 | { |
---|
| 1073 | return semicListMilnorWrong; |
---|
| 1074 | } |
---|
| 1075 | |
---|
| 1076 | // ------------------------- |
---|
| 1077 | // check geometrical genus |
---|
| 1078 | // ------------------------- |
---|
| 1079 | |
---|
| 1080 | for( pg=0, i=0; i<n; i++ ) |
---|
| 1081 | { |
---|
| 1082 | if( (*num)[i]<=(*den)[i] ) |
---|
| 1083 | { |
---|
| 1084 | pg += (*mul)[i]; |
---|
| 1085 | } |
---|
| 1086 | } |
---|
| 1087 | |
---|
| 1088 | if( pg != (int)(l->m[1].Data( )) ) |
---|
| 1089 | { |
---|
| 1090 | return semicListPGWrong; |
---|
| 1091 | } |
---|
| 1092 | |
---|
| 1093 | return semicOK; |
---|
| 1094 | } |
---|
| 1095 | |
---|
| 1096 | // ---------------------------------------------------------------------------- |
---|
| 1097 | // print out an error message for a spectrum list |
---|
| 1098 | // ---------------------------------------------------------------------------- |
---|
| 1099 | |
---|
| 1100 | void list_error( semicState state ) |
---|
| 1101 | { |
---|
| 1102 | switch( state ) |
---|
| 1103 | { |
---|
| 1104 | case semicListTooShort: |
---|
| 1105 | WerrorS( "the list is too short" ); |
---|
| 1106 | break; |
---|
| 1107 | case semicListTooLong: |
---|
| 1108 | WerrorS( "the list is too long" ); |
---|
| 1109 | break; |
---|
| 1110 | |
---|
| 1111 | case semicListFirstElementWrongType: |
---|
| 1112 | WerrorS( "first element of the list should be int" ); |
---|
| 1113 | break; |
---|
| 1114 | case semicListSecondElementWrongType: |
---|
| 1115 | WerrorS( "second element of the list should be int" ); |
---|
| 1116 | break; |
---|
| 1117 | case semicListThirdElementWrongType: |
---|
| 1118 | WerrorS( "third element of the list should be int" ); |
---|
| 1119 | break; |
---|
| 1120 | case semicListFourthElementWrongType: |
---|
| 1121 | WerrorS( "fourth element of the list should be intvec" ); |
---|
| 1122 | break; |
---|
| 1123 | case semicListFifthElementWrongType: |
---|
| 1124 | WerrorS( "fifth element of the list should be intvec" ); |
---|
| 1125 | break; |
---|
| 1126 | case semicListSixthElementWrongType: |
---|
| 1127 | WerrorS( "sixth element of the list should be intvec" ); |
---|
| 1128 | break; |
---|
| 1129 | |
---|
| 1130 | case semicListNNegative: |
---|
| 1131 | WerrorS( "first element of the list should be positive" ); |
---|
| 1132 | break; |
---|
| 1133 | case semicListWrongNumberOfNumerators: |
---|
| 1134 | WerrorS( "wrong number of numerators" ); |
---|
| 1135 | break; |
---|
| 1136 | case semicListWrongNumberOfDenominators: |
---|
| 1137 | WerrorS( "wrong number of denominators" ); |
---|
| 1138 | break; |
---|
| 1139 | case semicListWrongNumberOfMultiplicities: |
---|
| 1140 | WerrorS( "wrong number of multiplicities" ); |
---|
| 1141 | break; |
---|
| 1142 | |
---|
| 1143 | case semicListMuNegative: |
---|
| 1144 | WerrorS( "the Milnor number should be positive" ); |
---|
| 1145 | break; |
---|
| 1146 | case semicListPgNegative: |
---|
| 1147 | WerrorS( "the geometrical genus should be nonnegative" ); |
---|
| 1148 | break; |
---|
| 1149 | case semicListNumNegative: |
---|
| 1150 | WerrorS( "all numerators should be positive" ); |
---|
| 1151 | break; |
---|
| 1152 | case semicListDenNegative: |
---|
| 1153 | WerrorS( "all denominators should be positive" ); |
---|
| 1154 | break; |
---|
| 1155 | case semicListMulNegative: |
---|
| 1156 | WerrorS( "all multiplicities should be positive" ); |
---|
| 1157 | break; |
---|
| 1158 | |
---|
| 1159 | case semicListNotSymmetric: |
---|
| 1160 | WerrorS( "it is not symmetric" ); |
---|
| 1161 | break; |
---|
| 1162 | case semicListNotMonotonous: |
---|
| 1163 | WerrorS( "it is not monotonous" ); |
---|
| 1164 | break; |
---|
| 1165 | |
---|
| 1166 | case semicListMilnorWrong: |
---|
| 1167 | WerrorS( "the Milnor number is wrong" ); |
---|
| 1168 | break; |
---|
| 1169 | case semicListPGWrong: |
---|
| 1170 | WerrorS( "the geometrical genus is wrong" ); |
---|
| 1171 | break; |
---|
| 1172 | |
---|
| 1173 | default: |
---|
| 1174 | WerrorS( "unspecific error" ); |
---|
| 1175 | break; |
---|
| 1176 | } |
---|
| 1177 | } |
---|
| 1178 | |
---|
| 1179 | // ---------------------------------------------------------------------------- |
---|
| 1180 | // this procedure is called from the interpreter |
---|
| 1181 | // ---------------------------------------------------------------------------- |
---|
| 1182 | // first = list of spectrum numbers |
---|
| 1183 | // second = list of spectrum numbers |
---|
| 1184 | // result = sum of the two lists |
---|
| 1185 | // ---------------------------------------------------------------------------- |
---|
| 1186 | |
---|
| 1187 | BOOLEAN spaddProc( leftv result,leftv first,leftv second ) |
---|
| 1188 | { |
---|
| 1189 | semicState state; |
---|
| 1190 | |
---|
| 1191 | // ----------------- |
---|
| 1192 | // check arguments |
---|
| 1193 | // ----------------- |
---|
| 1194 | |
---|
| 1195 | lists l1 = (lists)first->Data( ); |
---|
| 1196 | lists l2 = (lists)second->Data( ); |
---|
| 1197 | |
---|
| 1198 | if( (state=list_is_spectrum( l1 )) != semicOK ) |
---|
| 1199 | { |
---|
| 1200 | WerrorS( "first argument is not a spectrum:" ); |
---|
| 1201 | list_error( state ); |
---|
| 1202 | } |
---|
| 1203 | else if( (state=list_is_spectrum( l2 )) != semicOK ) |
---|
| 1204 | { |
---|
| 1205 | WerrorS( "second argument is not a spectrum:" ); |
---|
| 1206 | list_error( state ); |
---|
| 1207 | } |
---|
| 1208 | else |
---|
| 1209 | { |
---|
| 1210 | spectrum s1( l1 ); |
---|
| 1211 | spectrum s2( l2 ); |
---|
| 1212 | spectrum sum( s1+s2 ); |
---|
| 1213 | |
---|
| 1214 | result->rtyp = LIST_CMD; |
---|
| 1215 | result->data = (char*)(sum.thelist( )); |
---|
| 1216 | } |
---|
| 1217 | |
---|
| 1218 | return (state!=semicOK); |
---|
| 1219 | } |
---|
| 1220 | |
---|
| 1221 | // ---------------------------------------------------------------------------- |
---|
| 1222 | // this procedure is called from the interpreter |
---|
| 1223 | // ---------------------------------------------------------------------------- |
---|
| 1224 | // first = list of spectrum numbers |
---|
| 1225 | // second = integer |
---|
| 1226 | // result = the multiple of the first list by the second factor |
---|
| 1227 | // ---------------------------------------------------------------------------- |
---|
| 1228 | |
---|
| 1229 | BOOLEAN spmulProc( leftv result,leftv first,leftv second ) |
---|
| 1230 | { |
---|
| 1231 | semicState state; |
---|
| 1232 | |
---|
| 1233 | // ----------------- |
---|
| 1234 | // check arguments |
---|
| 1235 | // ----------------- |
---|
| 1236 | |
---|
| 1237 | lists l = (lists)first->Data( ); |
---|
| 1238 | int k = (int)second->Data( ); |
---|
| 1239 | |
---|
| 1240 | if( (state=list_is_spectrum( l ))!=semicOK ) |
---|
| 1241 | { |
---|
| 1242 | WerrorS( "first argument is not a spectrum" ); |
---|
| 1243 | list_error( state ); |
---|
| 1244 | } |
---|
| 1245 | else if( k < 0 ) |
---|
| 1246 | { |
---|
| 1247 | WerrorS( "second argument should be positive" ); |
---|
| 1248 | state = semicMulNegative; |
---|
| 1249 | } |
---|
| 1250 | else |
---|
| 1251 | { |
---|
| 1252 | spectrum s( l ); |
---|
| 1253 | spectrum product( k*s ); |
---|
| 1254 | |
---|
| 1255 | result->rtyp = LIST_CMD; |
---|
| 1256 | result->data = (char*)product.thelist( ); |
---|
| 1257 | } |
---|
| 1258 | |
---|
| 1259 | return (state!=semicOK); |
---|
| 1260 | } |
---|
| 1261 | |
---|
| 1262 | // ---------------------------------------------------------------------------- |
---|
| 1263 | // this procedure is called from the interpreter |
---|
| 1264 | // ---------------------------------------------------------------------------- |
---|
| 1265 | // first = list of spectrum numbers |
---|
| 1266 | // second = list of spectrum numbers |
---|
| 1267 | // result = semicontinuity index |
---|
| 1268 | // ---------------------------------------------------------------------------- |
---|
| 1269 | |
---|
[a7f305] | 1270 | BOOLEAN semicProc3 ( leftv res,leftv u,leftv v,leftv w ) |
---|
[7885020] | 1271 | { |
---|
[5c2863] | 1272 | semicState state; |
---|
[cd6dbb2] | 1273 | BOOLEAN qh=(((int)w->Data())==1); |
---|
[7885020] | 1274 | |
---|
[5c2863] | 1275 | // ----------------- |
---|
| 1276 | // check arguments |
---|
| 1277 | // ----------------- |
---|
[7885020] | 1278 | |
---|
[cd6dbb2] | 1279 | lists l1 = (lists)u->Data( ); |
---|
| 1280 | lists l2 = (lists)v->Data( ); |
---|
[7885020] | 1281 | |
---|
[5c2863] | 1282 | if( (state=list_is_spectrum( l1 ))!=semicOK ) |
---|
| 1283 | { |
---|
| 1284 | WerrorS( "first argument is not a spectrum" ); |
---|
| 1285 | list_error( state ); |
---|
| 1286 | } |
---|
| 1287 | else if( (state=list_is_spectrum( l2 ))!=semicOK ) |
---|
| 1288 | { |
---|
| 1289 | WerrorS( "second argument is not a spectrum" ); |
---|
| 1290 | list_error( state ); |
---|
| 1291 | } |
---|
| 1292 | else |
---|
| 1293 | { |
---|
| 1294 | spectrum s1( l1 ); |
---|
| 1295 | spectrum s2( l2 ); |
---|
[7885020] | 1296 | |
---|
[cd6dbb2] | 1297 | res->rtyp = INT_CMD; |
---|
[5c2863] | 1298 | if (qh) |
---|
[cd6dbb2] | 1299 | res->data = (void*)(s1.mult_spectrumh( s2 )); |
---|
[5c2863] | 1300 | else |
---|
[cd6dbb2] | 1301 | res->data = (void*)(s1.mult_spectrum( s2 )); |
---|
[5c2863] | 1302 | } |
---|
[7885020] | 1303 | |
---|
[5c2863] | 1304 | // ----------------- |
---|
| 1305 | // check status |
---|
| 1306 | // ----------------- |
---|
[7885020] | 1307 | |
---|
[5c2863] | 1308 | return (state!=semicOK); |
---|
[7885020] | 1309 | } |
---|
[a7f305] | 1310 | BOOLEAN semicProc ( leftv res,leftv u,leftv v ) |
---|
[efaa52] | 1311 | { |
---|
[cd6dbb2] | 1312 | sleftv tmp; |
---|
| 1313 | memset(&tmp,0,sizeof(tmp)); |
---|
[a7f305] | 1314 | tmp.rtyp=INT_CMD; |
---|
| 1315 | /* tmp.data = (void *)0; -- done by memset */ |
---|
[efaa52] | 1316 | |
---|
[a7f305] | 1317 | return semicProc3(res,u,v,&tmp); |
---|
[efaa52] | 1318 | } |
---|
[7885020] | 1319 | #endif /* HAVE_SPECTRUM */ |
---|
| 1320 | // ---------------------------------------------------------------------------- |
---|
| 1321 | // spectrum.cc |
---|
| 1322 | // end of file |
---|
| 1323 | // ---------------------------------------------------------------------------- |
---|