Changeset 1d5418 in git
 Timestamp:
 Dec 15, 2004, 8:02:33 PM (20 years ago)
 Branches:
 (u'fiekerDuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
 Children:
 501f7c553e083fd3e3bfad340f57a827228a7911
 Parents:
 55e3087801d49cf32ed31572abf9343bd66a4d0a
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/center.lib
r55e308 r1d5418 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: center.lib,v 1.1 1 20041207 16:02:02plural Exp $";2 version="$Id: center.lib,v 1.12 20041215 19:02:33 plural Exp $"; 3 3 category="Noncommutative"; 4 4 info=" 5 5 LIBRARY: center.lib computation of central elements of Galgebras and Qalgebras. 6 AUTHOR S: Oleksandr Motsak, motsak@mathematik.unikl.de.6 AUTHOR: Oleksandr Motsak, motsak@mathematik.unikl.de. 7 7 OVERVIEW: 8 8 This is a library for computing the central elements and centralisators of elements in various noncommutative algebras. 9 9 Implementation is mainly based on algorithms, written in the frame of the diploma thesis by Oleksandr Motsak (advisor: Prof. S.A. Ovsienko, using the definitions from diploma thesis by V.Levandovskyy), at Kyiv Taras Shevchenko University: 'An algorithm for the computation of the center of noncommutative polynomial algebra'. 10 @* The last version of this library can be found via internet .10 @* The last version of this library can be found via internet at author's home page. 11 11 12 12 SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' … … 23 23 24 24 25 / //////////////////////////////////////////////////////////////////////////////25 /******************************************************/ 26 26 // stuff 27 27 28 / //////////////////////////////////////////////////////////////////////////////28 /******************************************************/ 29 29 static proc myValue ( def s, list # ) 30 "30 " 31 31 return s or (typeof(s))(#) 32 32 " 33 33 { 34 35 36 37 38 39 40 41 42 43 }; 44 45 / //////////////////////////////////////////////////////////////////////////////34 def @p = s; 35 if ( size(#) > 0 ) 36 { 37 if ( typeof( #[1] ) == typeof(s) ) 38 { 39 @p = #[1]; 40 }; 41 }; 42 return (@p); 43 }; 44 45 /******************************************************/ 46 46 static proc myInt ( list # ) 47 "47 " 48 48 return 0 or int(#) 49 49 " 50 50 { 51 52 53 }; 54 55 / //////////////////////////////////////////////////////////////////////////////51 int @p = 0; 52 return (myValue ( @p, # )); 53 }; 54 55 /******************************************************/ 56 56 static proc myPoly ( list # ) 57 "57 " 58 58 return 0 or poly(#) 59 59 " 60 60 { 61 62 63 }; 64 65 / //////////////////////////////////////////////////////////////////////////////61 poly @p = 0; 62 return (myValue( @p, # )); 63 }; 64 65 /******************************************************/ 66 66 static proc myRing ( list # ) 67 "67 " 68 68 return basring or ring(#) 69 69 " 70 70 { 71 72 73 }; 74 75 / //////////////////////////////////////////////////////////////////////////////71 def @p = basering; 72 return (myValue ( @p, # )); 73 }; 74 75 /******************************************************/ 76 76 static proc myString ( list # ) 77 "77 " 78 78 return basring or ring(#) 79 79 " 80 80 { 81 82 83 }; 84 85 / //////////////////////////////////////////////////////////////////////////////81 string @p = ""; 82 return (myValue ( @p, # )); 83 }; 84 85 /******************************************************/ 86 86 static proc myIdeal ( list # ) 87 "87 " 88 88 return 0 or ideal(#) 89 89 " 90 90 { 91 92 93 }; 94 95 / //////////////////////////////////////////////////////////////////////////////96 // /for debug97 98 / //////////////////////////////////////////////////////////////////////////////91 ideal @p; 92 return (myValue ( @p, # )); 93 }; 94 95 /******************************************************/ 96 // for debug 97 98 /******************************************************/ 99 99 static proc toprint() 100 100 { 101 102 }; 103 104 / //////////////////////////////////////////////////////////////////////////////105 static proc Print( def a)106 { 107 108 109 print (a);110 101 return (0); 102 }; 103 104 /******************************************************/ 105 static proc Print( list # ) 106 { 107 if ( toprint() ) 108 { 109 print (#); 110 }; 111 111 }; 112 112 113 / //////////////////////////////////////////////////////////////////////////////113 /******************************************************/ 114 114 static proc BCall(string Name, list #) 115 115 { 116 117 118 119 120 }; 121 122 / //////////////////////////////////////////////////////////////////////////////116 if( toprint() ) 117 { 118 "CALL: ", Name, "( ", string(#), " )"; 119 }; 120 }; 121 122 /******************************************************/ 123 123 static proc ECall(string Name, list #) 124 124 { 125 if( toprint() ) 126 { 127 "RET: ", Name, "(...)={ ", string(#), " }"; 128 }; 129 }; 130 131 132 //////////////////////////////////////////////////////////////////// 133 //////////////////////////////////////////////////////////////////// 125 if( toprint() ) 126 { 127 "RET : ", Name, "(...)={ ", string(#), " }"; 128 }; 129 }; 130 131 132 /******************************************************/ 134 133 // "my" degres functions and lie procedures... 135 134 136 / ///////////////////////////////////////////////////////////////////137 static proc maxDeg( proc givendeg,def z )138 "135 /******************************************************/ 136 static proc maxDeg( def z ) 137 " 139 138 returns: int : max of givendeg( z_i ), among all z_i \in z 140 139 returns 1 if z is empty or contain only 0 polynomial 141 140 " 142 141 { 143 int max = 1; int d; 144 145 for(int i=size(z); i>0; i ) 146 { 147 d = givendeg(z[i]); 148 if( d > max ) 149 { 150 max = d; 151 }; 152 }; 153 154 return(max); 155 }; 156 157 //////////////////////////////////////////////////////////////////// 158 static proc minDeg( proc givendeg, def z ) 159 " 160 returns: int : min of givendeg( z_i ), among all z_i \in z 161 returns 1 if z is empty or contain only 0 polynomial 162 " 163 { 164 if( size(z) == 0 ) 165 { 166 return(1); 167 }; 168 169 int min = 2; int d; 170 171 for(int i = size(z); i>0; i ) 172 { 173 d = givendeg(z[i]); 174 if( (min == 2) or (d < min) ) 175 { 176 min = d; 177 }; 178 179 }; 180 181 return(min); 182 }; 183 184 //////////////////////////////////////////////////////////////////// 185 static proc myMonomDeg( poly p) 186 " 187 input: poly p 188 return: just number of veriables in lead monomial of p! 189 Note: for use by maxList and minList 190 " 191 { 192 if( p == 0 ) 193 { 194 return (1); 195 }; 196 197 intvec e = leadexp(p); 198 int d = 0; 199 200 for ( int i = size(e); i>0; i ) 201 { 202 d = d + e[i]; 203 }; 204 205 return (d); 206 }; 207 208 //////////////////////////////////////////////////////////////////// 209 static proc stdDeg ( poly p) 210 " 211 input: poly p 212 return: degree of p 213 Note: just in order to have it as 'proc' (for use by maxList) 214 " 215 { 216 return (deg(p)); 217 }; 218 219 //////////////////////////////////////////////////////////////////// 220 static proc myDeg ( poly p) 221 " 222 input: poly p 223 return: 'my degree' of p 224 no real use... just for testing... 225 " 226 { 227 return( maxDeg(myMonomDeg,p) ); 228 }; 229 230 /////////////////////////////////////////////////////////////////////////////// 231 /////////////////////////////////////////////////////////////////////////////// 232 // some bisiness for my_var and lie_rank 233 234 //////////////////////////////////////////////////////////////////// 142 int max = 1; 143 int d; 144 145 for(int i=size(z); i>0; i ) 146 { 147 d = deg(z[i]); 148 if( d > max ) 149 { 150 max = d; 151 }; 152 }; 153 154 return(max); 155 }; 156 157 /******************************************************/ 158 // some bisiness for my_var 159 160 /******************************************************/ 235 161 static proc myCoeff ( poly f, poly m ) 236 "162 " 237 163 input: poly f, 238 164 return: coeff at m 239 165 " 240 166 { 241 m = leadmonom(m); 242 243 int i = size(f); 244 245 while ( (i>0) and (leadmonom(f[i])<m) ) 246 { 247 i; 248 }; 249 if( i == 0 ) 250 { 251 return ( 0 ); 252 }; 253 254 if(leadmonom(f[i]) == m) 255 { 256 return ( leadcoef(f[i]) ); 257 }; 258 259 return ( 0 ); 260 }; 261 262 263 /////////////////////////////////////////////////////////////////////////////// 167 m = leadmonom(m); 168 169 int i = size(f); 170 171 while ( (i>0) and (leadmonom(f[i])<m) ) 172 { 173 i; 174 }; 175 if( i == 0 ) 176 { 177 return ( 0 ); 178 }; 179 180 if(leadmonom(f[i]) == m) 181 { 182 return ( leadcoef(f[i]) ); 183 }; 184 185 return ( 0 ); 186 }; 187 188 /******************************************************/ 264 189 static proc my_vars () 265 190 { 266 267 268 269 270 271 272 273 274 275 }; 276 277 / //////////////////////////////////////////////////////////////////////////////191 int N = nvars( basering ); 192 list V = list(); 193 194 for ( int k = 1; k <= N; k++ ) 195 { 196 V[k] = var(k); 197 }; 198 199 return (V); 200 }; 201 202 /******************************************************/ 278 203 static proc my_commutators () 279 204 { 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 }; 300 301 / //////////////////////////////////////////////////////////////////////////////205 list V = my_vars(); 206 int N = size( V ); 207 208 matrix M[N][N] = 0; 209 210 poly v, d; int i; 211 212 for ( int k = 2; k <= N; k++ ) 213 { 214 v = V[k]; 215 for ( i = 1; i < k; i++ ) 216 { 217 d = V[i]; 218 M[k,i] = v*d  d*v; // [var(k),var(i)] 219 M[i,k] = M[k,i]; // [var(i),var(k)] == [var(k),var(i)] 220 }; 221 }; 222 223 return (M); 224 }; 225 226 /******************************************************/ 302 227 static proc my_associated () 303 228 " … … 308 233 " 309 234 { 310 int i, j; poly a; poly d; 311 312 int N = nvars( basering ); 313 list V = my_vars(); 314 matrix M = my_commutators(); 315 316 matrix A[N][N]; 317 318 int cartan; 319 320 list RESULT = list(); 321 int r_begin = 1; 322 int r_end = N; 323 324 325 for ( int k = 1; k <= N; k++ ) 326 { 327 A = 0; 328 cartan = 1; 329 for ( i = 1; i <= N; i++ ) 235 int i, j; poly a; poly d; 236 237 int N = nvars( basering ); 238 list V = my_vars(); 239 matrix M = my_commutators(); 240 241 matrix A[N][N]; 242 243 int cartan; 244 245 list RESULT = list(); 246 int r_begin = 1; 247 int r_end = N; 248 249 250 for ( int k = 1; k <= N; k++ ) 251 { 252 A = 0; 253 cartan = 1; 254 for ( i = 1; i <= N; i++ ) 255 { 256 d = M[k,i]; 257 for ( j = 1; j <= N; j++ ) 258 { 259 a = myCoeff( d, V[j] ); 260 A[i,j] = a; 261 262 if( (i!=j) and (a!=0) ) 330 263 { 331 d = M[k,i]; 332 for ( j = 1; j <= N; j++ ) 333 { 334 a = myCoeff( d, V[j] ); A[i,j] = a; 335 336 if( (i!=j) and (a!=0) ) 337 { 338 cartan = 0; 339 }; 340 }; 341 }; 342 343 if ( cartan ) 344 { 345 RESULT[r_begin] = list( V[k], cartan, A ); 346 r_begin++; 347 } else 348 { 349 RESULT[r_end] = list( V[k], cartan, A ); 350 r_end; 351 }; 264 cartan = 0; 265 }; 266 }; 267 }; 268 269 if ( cartan ) 270 { 271 RESULT[r_begin] = list( V[k], cartan, A ); 272 r_begin++; 273 } else 274 { 275 RESULT[r_end] = list( V[k], cartan, A ); 276 r_end; 277 }; 352 278 353 354 355 356 357 }; 358 359 /****************************************************** ************************/279 }; 280 281 return (RESULT); 282 283 }; 284 285 /******************************************************/ 360 286 static proc my_var_init() 361 287 " … … 364 290 { 365 291 if( defined( @@@SORTEDVARARRAY ) ) 366 {367 368 };292 { 293 kill ( @@@SORTEDVARARRAY ); 294 }; 369 295 370 296 list @@@SORTEDVARARRAY; … … 374 300 375 301 for ( int i = 1; i<= size(V); i++ ) 376 {302 { 377 303 @@@SORTEDVARARRAY[i] = V[i][1]; 378 };304 }; 379 305 380 306 Print( "@@@SORTEDVARARRAY: " + string(@@@SORTEDVARARRAY) ); … … 383 309 }; 384 310 385 /****************************************************** ************************/311 /******************************************************/ 386 312 static proc my_var(int @number ) 387 "313 " 388 314 'my' var (with my order of variables), 389 315 before use call my_var_init() … … 391 317 { 392 318 if( defined( @@@SORTEDVARARRAY ) ) 393 {394 395 };319 { 320 return( @@@SORTEDVARARRAY[@number] ); 321 }; 396 322 397 323 Print( "Error: my_var_init() was not called before this..." ); … … 400 326 }; 401 327 402 /****************************************************** ************************/328 /******************************************************/ 403 329 static proc my_var_done() 404 "330 " 405 331 this procedure should be called after finishing using my_var in order to clean up. 406 332 " 407 333 { 408 334 if( defined( @@@SORTEDVARARRAY ) ) 409 { 410 kill ( @@@SORTEDVARARRAY ); 411 }; 412 }; 413 414 415 416 /////////////////////////////////////////////////////////////////////////////// 417 static proc index( int i ) 418 " 419 i need indexing of arrays from 0 => ARRAY[ index(i) ]  everywhere... 420 " 421 { 422 return (i + 1); 423 }; 424 425 /////////////////////////////////////////////////////////////////////////////// 335 { 336 kill ( @@@SORTEDVARARRAY ); 337 }; 338 }; 339 340 341 /******************************************************/ 342 // QNF stuff 343 344 /******************************************************/ 426 345 static proc myBitParam( int param, int n ) 427 "346 " 428 347 return: string: param in bin within no less then n digits (adding zeroes) 429 348 " 430 349 { 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 }; 446 447 / //////////////////////////////////////////////////////////////////////////////350 string s = ""; 351 352 while ( param != 0 ) 353 { 354 s = string ( param % 2 ) + s; 355 param = param / 2; 356 n ; 357 }; 358 while ( n > 0 ) 359 { 360 s = "0" + s; 361 n ; 362 }; 363 return ( s ); 364 }; 365 366 /******************************************************/ 448 367 static proc QNF_init(list #) 449 368 { 450 451 { 452 453 }; 454 455 456 { 457 458 /// setup redSB, redTail options459 460 461 462 463 464 }; 465 }; 466 467 / //////////////////////////////////////////////////////////////////////////////369 if( defined(@@@MY_QIDEAL) ) 370 { 371 kill (@@@MY_QIDEAL); 372 }; 373 374 if( typeof(basering) == "qring" ) 375 { 376 ideal @@@MY_QIDEAL = twostd(myIdeal(#)); 377 // setup redSB, redTail options 378 export(@@@MY_QIDEAL); 379 380 option(redSB); 381 option(redTail); 382 383 }; 384 }; 385 386 /******************************************************/ 468 387 static proc QNF_done() 469 388 { 470 471 { 472 473 }; 474 }; 475 476 / //////////////////////////////////////////////////////////////////////////////477 static proc QNF_poly ( poly p, list # ) )478 { 479 480 { 481 /*389 if( defined(@@@MY_QIDEAL) ) 390 { 391 kill (@@@MY_QIDEAL); 392 }; 393 }; 394 395 /******************************************************/ 396 static proc QNF_poly ( poly p, list # ) 397 { 398 if( defined(@@@MY_QIDEAL) ) 399 { 400 /* 482 401 int param = 1 ; // by def: reduce both 1st and 2nd entries 483 402 param = myValue(param, #); 484 403 string bits = myBitParam( param, 8 ); 485 */486 487 488 489 490 491 }; // /QRing492 493 494 }; 495 496 497 / //////////////////////////////////////////////////////////////////////////////404 */ 405 BCall( "QNF_poly", p ); 406 407 p = NF( p, @@@MY_QIDEAL ); 408 409 ECall( "QNF_poly", p ); 410 }; // QRing 411 412 return(p); 413 }; 414 415 416 /******************************************************/ 498 417 static proc QNF_list ( list l, list # ) 499 "418 " 500 419 expects list of records with two polynomial entries 501 420 and in case of Qring allpies NF( * , std(0) ) to entries (depends on parameter) … … 516 435 { 517 436 518 if( defined(@@@MY_QIDEAL) ) 519 { 520 521 int param = 1 + 2; // by def: reduce both 1st and 2nd entries 522 param = myValue(param, #); 523 524 string bits = myBitParam( param, 8 ); 525 526 BCall( "QNF_list", "list, " + bits ); 527 528 poly temp; 529 530 for ( int i = size(l); i>0 ; i  ) 531 { 532 533 if ( typeof( l[i] ) == "poly" ) 534 { 535 536 if ( (bits[8] == "1") or (bits[6] == "1") ) 437 if( defined(@@@MY_QIDEAL) ) 438 { 439 440 int param = 1 + 2; // by def: reduce both 1st and 2nd entries 441 param = myValue(param, #); 442 443 string bits = myBitParam( param, 8 ); 444 445 BCall( "QNF_list", "list, " + bits ); 446 447 poly temp; 448 449 for ( int i = size(l); i>0 ; i  ) 450 { 451 452 if ( typeof( l[i] ) == "poly" ) 453 { 454 455 if ( (bits[8] == "1") or (bits[6] == "1") ) 456 { 457 temp = NF( l[i], @@@MY_QIDEAL ); 458 459 if ( bits[6] == "1" ) 460 {// for PBW in qrings: kill out pbw entries where pbw monom was reduced 461 462 if( temp != l[i] ) 463 { 464 l = delete( l, i ); 465 i ; 466 continue; 467 }; 468 }; 469 470 if ( bits[8] == "1" ) 537 471 { 538 temp = NF( l[i], @@@MY_QIDEAL ); 472 l[i] = temp; 473 }; 474 }; 475 }; 476 477 if ( typeof( l[i] ) == "list" ) 478 { 479 // 1st 480 481 if ( size(l[i])>0 ) 482 { 483 if ( (bits[8] == "1") or (bits[6] == "1") ) 484 { 485 if( typeof( l[i][1] ) == "poly" ) 486 { 539 487 540 if ( bits[6] == "1" ) 541 {/// for PBW in qrings: kill out pbw entries where pbw monom was reduced 488 temp = NF( l[i][1], @@@MY_QIDEAL ); 542 489 543 if( temp != l[i] ) 490 if ( bits[6] == "1" ) 491 {// for PBW in qrings: kill out pbw entries where pbw monom was reduced 492 493 if( temp != l[i][1] ) 544 494 { 545 546 547 495 l = delete( l, i ); 496 i ; 497 continue; 548 498 }; 549 499 }; 550 500 551 501 if ( bits[8] == "1" ) 552 502 { 553 l[i] = temp; 554 }; 555 }; 556 }; 557 558 if ( typeof( l[i] ) == "list" ) 559 { 560 /// 1st 561 562 if ( size(l[i])>0 ) 563 { 564 if ( (bits[8] == "1") or (bits[6] == "1") ) 565 { 566 if( typeof( l[i][1] ) == "poly" ) 567 { 568 569 temp = NF( l[i][1], @@@MY_QIDEAL ); 570 571 if ( bits[6] == "1" ) 572 {/// for PBW in qrings: kill out pbw entries where pbw monom was reduced 573 574 if( temp != l[i][1] ) 575 { 576 l = delete( l, i ); 577 i ; 578 continue; 579 }; 503 l[i][1] = temp; 580 504 }; 581 505 582 if ( bits[8] == "1" )583 {584 l[i][1] = temp;585 };586 587 506 }; 588 507 }; 589 508 }; 590 591 509 510 // 2nd 592 511 593 512 if ( size(l[i])>1 ) 594 513 { 595 514 if ( bits[7] == "1" ) 596 515 { 597 516 if( typeof( l[i][2] ) == "poly" ) 598 517 { 599 600 601 518 temp = NF( l[i][2], @@@MY_QIDEAL ); 519 520 l[i][2] = temp; 602 521 }; 603 522 }; … … 605 524 }; 606 525 }; 607 608 609 526 527 ECall( "QNF_list", "list" ); 528 610 529 }; // Qring 611 530 612 return ( l ); 613 }; 614 615 /////////////////////////////////////////////////////////////////////////////// 616 /// things for zCommon 617 618 /******************************************************************************/ 531 return ( l ); 532 }; 533 534 /******************************************************/ 535 // things for zCommon 536 537 /******************************************************/ 538 static proc index( int i ) 539 " 540 i need indexing of arrays from 0 => ARRAY[ index(i) ]  everywhere... 541 " 542 { 543 return (i + 1); 544 }; 545 546 /******************************************************/ 619 547 static proc uni_poly( poly p ) 620 "548 " 621 549 returns polynomial with the same monomials but without coefficients. 622 550 " 623 551 { 624 625 626 627 628 629 630 }; 631 632 /****************************************************** ************************/552 poly @tt = poly(0); 553 for ( int @k = size(p); @k > 0; @k  ) 554 { 555 @tt = @tt + leadmonom(p[@k]); 556 }; 557 return (@tt); 558 }; 559 560 /******************************************************/ 633 561 static proc get_max ( poly @t, number @def ) 634 562 { 635 int @n = size( @t ); 563 int @n = size( @t ); 564 565 if ( @n == 0 ) 566 { 567 return (@def); 568 }; 636 569 637 if ( @n == 0 ) 638 { 639 return (@def); 640 }; 641 642 number @max = leadcoef(@t[1]); 643 number @mm; 644 645 if ( @n > 1) 646 { 647 for ( int @i = 2; @i <= @n ;@i ++ ) 648 { 649 @mm = leadcoef ( @t[@i] ); 650 if ( @mm < 0 ) 651 { 652 @mm = @mm; 653 }; 654 655 if( @mm > @max ) 656 { 657 @max = @mm; 658 }; 659 }; 660 }; 661 662 @max = @max + 1; 663 if ( @max == 0 ) 664 { 665 @max = @max + 1; 666 }; 667 return ( 1/@max ); 668 }; 669 670 671 /******************************************************************************/ 570 number @max = leadcoef(@t[1]); 571 number @mm; 572 573 if ( @n > 1) 574 { 575 for ( int @i = 2; @i <= @n ;@i ++ ) 576 { 577 @mm = leadcoef ( @t[@i] ); 578 if ( @mm < 0 ) 579 { 580 @mm = @mm; 581 }; 582 583 if( @mm > @max ) 584 { 585 @max = @mm; 586 }; 587 }; 588 }; 589 590 @max = @max + 1; 591 if ( @max == 0 ) 592 { 593 @max = @max + 1; 594 }; 595 return ( 1/@max ); 596 }; 597 598 599 /******************************************************/ 672 600 static proc get_uni_sum(list @given, int @num) 673 601 { 674 675 676 677 678 679 680 681 682 { 683 684 { 685 602 int @l, @k; 603 604 poly @t, @tt; 605 606 @l = size (@given); 607 608 @t = poly(0); 609 for ( @k = @l; @k > 0; @k  ) 610 { 611 if (@num == 1) 612 { 613 @tt = @given[@k]; 686 614 } else 687 {615 { 688 616 @tt = @given[@k][2]; 689 };690 691 692 }; 693 694 695 }; 696 697 698 699 / ///////////////////////////////////////////////////////////////////617 }; 618 619 @t = @t + uni_poly( @tt ); 620 }; 621 622 return ( uni_poly(@t) ); 623 }; 624 625 626 627 /******************************************************/ 700 628 static proc LM ( intvec exp ) 701 "629 " 702 630 input : given exponent 703 631 return: monom with this exponent... 704 632 " 705 633 { 706 poly @f = 1; int @deg; 707 for( int @i = 1; @i <= size(exp); @i ++ ) 708 { 709 @deg = exp[@i]; 710 if ( @deg > 0 ) 711 { 712 @f = @f * (var(@i)^(@deg)); 713 }; 714 }; 715 716 return (@f); 717 }; 718 719 720 721 722 ///////////////////////////////////////////////////////////////////////////// 723 ///////////////////////////////////////////////////////////////////////////// 724 /// reduced centers  computation of "minimal" set of generators 725 726 LIB "general.lib"; /// for sort 727 728 ////////////////////////////////////////////////////////////////////////// 634 poly @f = 1; int @deg; 635 for( int @i = 1; @i <= size(exp); @i ++ ) 636 { 637 @deg = exp[@i]; 638 if ( @deg > 0 ) 639 { 640 @f = @f * (var(@i)^(@deg)); 641 }; 642 }; 643 644 return (@f); 645 }; 646 647 /******************************************************/ 648 // reduced centers  computation of "minimal" set of generators 649 650 LIB "general.lib"; // for sort 651 652 /******************************************************/ 729 653 static proc zSort ( list @z ) 730 "654 " 731 655 Sort elements of a list of polynoms, 732 656 and normalize 733 657 " 734 658 { 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 { /// all zeroes766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 p = p* (1/leadcoef(p)); /// normalize796 797 798 799 // "OUT SORT::";800 // result;801 802 803 }; 804 805 806 807 / /////////////////////////////////////////////////////////////////////////659 int n = size(@z); 660 if ( n < 1 ) 661 {// empty list 662 return( @z ); 663 }; 664 665 if ( n == 1 ) 666 { 667 if ( @z[1] == 0 ) // if zero => empty list 668 { 669 return(list()); 670 }; 671 672 @z[1] = @z[1] * ( 1/leadcoef(@z[1]) ) ; 673 674 return( @z ); 675 }; 676 677 int i = 1; 678 679 while ( i<=n ) 680 { 681 if (size( @z[i] ) != 0) 682 { 683 break; 684 }; 685 i++; 686 }; 687 688 if ( i > n ) 689 { // all zeroes 690 return(list()); 691 }; 692 693 ideal id = @z[i]; 694 i++; 695 696 while ( i<= n ) 697 { 698 if( @z[i] != 0 ) 699 { 700 id=@z[i],id; 701 }; 702 i++; 703 }; 704 705 // can sort only ideal generators:( 706 list srt = sort(id); 707 708 poly p; 709 // convert srt[1]:ideal to result:list. 710 list result = list(); 711 for ( i = size(srt[1]); i>0; i ) 712 { 713 p = srt[1][i]; 714 if ( p == 0 ) 715 { 716 i ; 717 continue; 718 }; 719 p = p* (1/leadcoef(p)); // normalize 720 result = list(p) + result; 721 }; 722 723 // "OUT SORT::"; 724 // result; 725 726 return ( result ); 727 }; 728 729 730 731 /******************************************************/ 808 732 static proc zRefine ( list @z ) 809 "733 " 810 734 kill terms by leading monomials... 811 735 Note: based on myCoeff 812 736 " 813 737 { 814 BCall("zRefine", "list"); 815 816 int @i, @j; 817 818 poly @f = 0; // leadmonom 819 820 poly @ff = 0; // polyes 821 poly @gg = 0; 822 number @nf, @ng; 823 824 int flag = 1; 825 826 while ( flag == 1 ) 827 { 828 flag = 0; 829 830 @z = zSort(@z); // sort out, < ... 831 832 if( size(@z) < 2 ) 738 BCall("zRefine", "list"); 739 740 int @i, @j; 741 742 poly @f = 0; // leadmonom 743 744 poly @ff = 0; // polyes 745 poly @gg = 0; 746 number @nf, @ng; 747 748 int flag = 1; 749 750 while ( flag == 1 ) 751 { 752 flag = 0; 753 754 @z = zSort(@z); // sort out, < ... 755 756 if( size(@z) < 2 ) 757 { 758 return (@z); 759 }; 760 761 for ( @i = size(@z); @i > 0; @i  ) // at 1st the biggest ... then smallest... 762 { 763 764 @ff = @z[@i]; 765 766 if( size(@ff) == 0 ) 767 { 768 @z = delete ( @z , @i ); 769 @i ; 770 continue; 771 }; 772 773 @ff = @ff*(1/leadcoef(@ff)); 774 @z[@i] = @ff; 775 @f = leadmonom(@ff); 776 777 for ( @j = (@i1); (@j>0); @j  ) 778 { 779 @gg = @z[@j]; 780 @ng = myCoeff(@gg, @f); // leads? 781 if( @ng!=0 ) 833 782 { 834 return (@z); 835 }; 836 837 for ( @i = size(@z); @i > 0; @i  ) // at 1st the biggest ... then smallest... 783 @z[@j] = @gg  @ng * @ff; 784 flag = 1; 785 }; 786 }; 787 for ( @j = (@i+1); (@j<=size(@z)); @j ++ ) 788 { 789 @gg = @z[@j]; 790 @ng = myCoeff(@gg, @f); 791 if( @ng!=0 ) 838 792 { 839 840 @ff = @z[@i]; 841 842 if( size(@ff) == 0 ) 843 { 844 @z = delete ( @z , @i ); 845 @i ; 846 continue; 847 }; 848 849 @ff = @ff*(1/leadcoef(@ff)); 850 @z[@i] = @ff; 851 @f = leadmonom(@ff); 852 853 for ( @j = (@i1); (@j>0); @j  ) 854 { 855 @gg = @z[@j]; 856 @ng = myCoeff(@gg, @f); /// leads? 857 if( @ng!=0 ) 858 { 859 @z[@j] = @gg  @ng * @ff; 860 flag = 1; 861 }; 862 }; 863 for ( @j = (@i+1); (@j<=size(@z)); @j ++ ) 864 { 865 @gg = @z[@j]; 866 @ng = myCoeff(@gg, @f); 867 if( @ng!=0 ) 868 { 869 @z[@j] = @gg  @ng * @ff; 870 flag = 1; 871 }; 872 }; 873 874 }; 875 }; 876 877 ECall("zRefine", "list"); 878 879 return ( @z ); 880 881 }; 882 883 884 //////////////////////////////////////////////////////////////////////////////// 885 /// procedures for building "bad leadmonomials" set 886 887 888 //////////////////////////////////////////////////////////////////////////////// 793 @z[@j] = @gg  @ng * @ff; 794 flag = 1; 795 }; 796 }; 797 798 }; 799 }; 800 801 ECall("zRefine", "list"); 802 803 return ( @z ); 804 805 }; 806 807 808 /******************************************************/ 809 // procedures for building "bad leadmonomials" set 810 811 812 /******************************************************/ 889 813 static proc checkPolyUniq( list l, poly p ) 890 "814 " 891 815 check wheather p sits already in l, assume l to be sizesorted < 892 816 return: 1 if present … … 894 818 " 895 819 { 896 897 898 899 900 901 902 903 904 905 906 907 820 BCall( "checkPolyUniq", "{ " + string(l) + " }, " + string(p) ); 821 int n = size(l); 822 int i = 1; 823 824 int s = size(p); 825 826 while( i<= n ) 827 { 828 if ( size(l[i]) >= s ) 829 { 830 break; 831 }; 908 832 909 910 911 912 /// now: size(l[i]) >= s913 914 915 916 917 833 i ++; 834 }; 835 836 // now: size(l[i]) >= s 837 while( i<= n ) 838 { 839 if ( size(l[i]) == s ) 840 { 841 break; 918 842 919 920 921 922 923 924 925 926 927 928 929 930 931 }; 932 933 / ///////////////////////////////////////////////////////////////////////////////843 }; 844 if ( l[i] == p ) 845 { 846 ECall( "checkPolyUniq", 1 ); 847 return (1); 848 }; 849 i ++; 850 }; 851 852 ECall( "checkPolyUniq", 1 ); 853 return ( 1 ); 854 855 }; 856 857 /******************************************************/ 934 858 static proc addPolyUniq( list l, poly p ) 935 "859 " 936 860 add p into l uniquely, and keep l sizesorted < 937 861 " 938 862 { 939 940 941 942 943 944 945 863 BCall( "addPolyUniq", " { " + string(l) + " }, " + string(p) ); 864 865 int n = size(l); 866 867 if ( n == 0 ) 868 { 869 l = list(p); 946 870 947 871 ECall( "addPolyUniq", l ); 948 872 949 950 951 952 953 954 955 956 957 958 959 960 961 873 return (l); 874 }; 875 876 int s = size(p); 877 878 int i = 1; 879 while( i<= n ) 880 { 881 if( size(l[i]) > s ) 882 { 883 l = insert( l, p, i1 ) ; 884 break; 885 }; 962 886 963 964 965 966 967 968 969 887 if( size(l[i]) == s ) 888 { 889 if ( l[i] == p ) 890 { 891 break; 892 }; 893 }; 970 894 971 972 973 974 975 976 977 978 979 980 981 }; 982 983 984 / ///////////////////////////////////////////////////////////////////////////////895 i++; 896 }; 897 898 if( i > n ) 899 { 900 l = l + list(p); 901 }; 902 903 ECall( "addPolyUniq", l ); 904 return(l); 905 }; 906 907 908 /******************************************************/ 985 909 static proc mergePolysUniq( list a, list b ) 986 "910 " 987 911 merge lists uniq 988 912 " 989 913 { 990 991 992 993 994 995 996 997 998 999 1000 }; 1001 1002 1003 / ///////////////////////////////////////////////////////////////////////////////914 BCall( "mergePolysUniq", "{ " + string(a) + " }, { " + string(b) + "} " ); 915 916 for( int i = 1; i <= size(b); i ++ ) 917 { 918 a = addPolyUniq(a, b[i]); 919 }; 920 921 ECall( "mergePolysUniq", a ); 922 923 return (a); 924 }; 925 926 927 /******************************************************/ 1004 928 static proc sortPolysUniq( list a ) 1005 "929 " 1006 930 sort list uniq 1007 931 " 1008 932 { 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 }; 1028 1029 / ///////////////////////////////////////////////////////////////////////////////933 BCall( "sortPolysUniq", a ); 934 935 if( size(a) < 2 ) 936 { 937 ECall( "sortPolysUniq", a ); 938 return(a); 939 }; 940 941 list b = list(a[1]); 942 943 for( int i = 2; i <= size(a); i ++ ) 944 { 945 b = addPolyUniq(b, a[i]); 946 }; 947 948 ECall( "sortPolysUniq", b ); 949 950 return (b); 951 }; 952 953 /******************************************************/ 1030 954 static proc addRecordUniq ( list leadD, list newD, intvec texp, poly tm, int kind, list prs ) 1031 "955 " 1032 956 if kind = 0 => for PBW  no products 1033 957 if kind = 1 => with products 1034 958 " 1035 959 { 1036 BCall ("addRecordUniq", "{old leads}, {new lead}, [" + string(texp) + "], " + string(tm) + ", " + string(kind) + ", {" + string(prs) + "}"); 1037 1038 int i; 1039 int f_add = 0; 1040 1041 prs = sortPolysUniq(prs); 1042 1043 /// trick: 1044 /// check for presens of a monomial @tm in current index poly of @leads (=> in list @leads) 1045 /// if char = 0 then new_size > (if not present) or = (if already present) 1046 /// if char > 0 then new_size > (if not present) or =< (if already present) 1047 /// !!!!! 1048 if( size(tm + leadD[2]) > size(leadD[2]) ) 1049 { 1050 f_add = 1; 960 BCall ("addRecordUniq", "{old leads}, {new lead}, [" + string(texp) + "], " + string(tm) + ", " + string(kind) + ", {" + string(prs) + "}"); 961 962 int i; 963 int f_add = 0; 964 965 prs = sortPolysUniq(prs); 966 967 // trick: 968 // check for presens of a monomial @tm in current index poly of @leads (=> in list @leads) 969 // if char = 0 then new_size > (if not present) or = (if already present) 970 // if char > 0 then new_size > (if not present) or =< (if already present) 971 // !!!!! 972 if( size(tm + leadD[2]) > size(leadD[2]) ) 973 { 974 f_add = 1; 975 } else 976 { 977 if ( kind != 0 ) 978 { 979 for ( i = 1; i<= size(leadD[1]); i++ ) 980 { 981 if ( leadD[1][i][2] == tm ) 982 { 983 for ( i = size(prs); i>0; i ) 984 { 985 f_add = checkPolyUniq( leadD[1][i][3], prs[i] ); 986 if( f_add == 1 ) 987 { 988 prs = delete(prs, i); 989 }; 990 }; 991 992 break; 993 }; 994 }; 995 }; 996 }; 997 998 if ( f_add == 1 ) 999 { 1000 1001 list newlist ; 1002 if ( kind != 0 ) 1003 { 1004 newlist = list ( list ( texp, tm, prs ) ); 1051 1005 } else 1052 { 1053 if ( kind != 0 ) 1006 { 1007 newlist = list ( list ( texp, tm ) ); 1008 }; 1009 1010 1011 if ( size(newD[1]) == 0 ) 1012 { 1013 newD[1] = newlist; 1014 newD[2] = tm; 1015 } else 1016 { 1017 1018 if( size(tm + newD[2]) > size(newD[2]) ) 1019 { 1020 newD[1] = newD[1] + newlist; 1021 newD[2] = newD[2] + tm; 1022 } else 1054 1023 { 1055 for ( i = 1; i<= size(leadD[1]); i++ ) 1024 if ( kind != 0 ) 1025 { 1026 for ( i = 1; i<= size(newD[1]); i++ ) 1056 1027 { 1057 if ( leadD[1][i][2] == tm ) 1028 if ( newD[1][i][2] == tm ) 1029 { 1030 newD[1][i][3] = mergePolysUniq( newD[1][i][3], prs ); 1031 break; 1032 }; 1033 }; 1034 }; 1035 }; 1036 }; 1037 }; 1038 1039 ECall("addRecordUniq", "{new leads}"); 1040 return (newD); 1041 }; 1042 1043 1044 /******************************************************/ 1045 static proc mergeRecordsUniq ( list old, list new, int kind ) 1046 { 1047 BCall ("mergeRecordsUniq", "{old leads}, {new leads}, " + string(kind) ); 1048 1049 if( size(new[1]) > 0 ) 1050 { 1051 if ( size (old[1]) == 0 ) 1052 { 1053 old = new; 1054 } else 1055 { 1056 if ( kind != 0 ) 1057 { 1058 int io; 1059 for ( int in = 1; in <= size( new[1] ); in ++ ) 1060 { 1061 if( size( new[1][in][2] + old[2] ) > size( old[2] ) ) 1062 { 1063 old[1] = old[1] + list(new[1][in]); 1064 old[2] = old[2] + new[1][in][2]; 1065 } else 1066 { 1067 for( io = 1; io <= size( old[1] ); io ++ ) 1068 { 1069 if( old[1][io][2] == new[1][in][2] ) 1058 1070 { 1059 for ( i = size(prs); i>0; i ) 1060 { 1061 f_add = checkPolyUniq( leadD[1][i][3], prs[i] ); 1062 if( f_add == 1 ) 1063 { 1064 prs = delete(prs, i); 1065 }; 1066 }; 1067 1068 break; 1071 old[1][io][3] = mergePolysUniq( old[1][io][3], new[1][in][3] ); 1072 break; 1069 1073 }; 1070 }; 1074 }; 1075 }; 1076 }; 1077 } else 1078 { 1079 old[1] = old[1] + new[1]; 1080 old[2] = old[2] + new[2]; 1071 1081 }; 1072 }; 1073 1074 if ( f_add == 1 ) 1075 { 1076 1077 list newlist ; 1078 if ( kind != 0 ) 1079 { 1080 newlist = list ( list ( texp, tm, prs ) ); 1081 } else 1082 { 1083 newlist = list ( list ( texp, tm ) ); 1084 }; 1085 1086 1087 if ( size(newD[1]) == 0 ) 1088 { 1089 newD[1] = newlist; 1090 newD[2] = tm; 1091 } else 1092 { 1093 1094 if( size(tm + newD[2]) > size(newD[2]) ) 1095 { 1096 newD[1] = newD[1] + newlist; 1097 newD[2] = newD[2] + tm; 1098 } else 1099 { 1100 if ( kind != 0 ) 1101 { 1102 for ( i = 1; i<= size(newD[1]); i++ ) 1103 { 1104 if ( newD[1][i][2] == tm ) 1105 { 1106 newD[1][i][3] = mergePolysUniq( newD[1][i][3], prs ); 1107 break; 1108 }; 1109 }; 1110 }; 1111 }; 1112 }; 1113 }; 1114 1115 ECall("addRecordUniq", "{new leads}"); 1116 return (newD); 1117 }; 1118 1119 1120 //////////////////////////////////////////////////////////////////////////////// 1121 static proc mergeRecordsUniq ( list old, list new, int kind ) 1122 { 1123 BCall ("mergeRecordsUniq", "{old leads}, {new leads}, " + string(kind) ); 1124 1125 if( size(new[1]) > 0 ) 1126 { 1127 if ( size (old[1]) == 0 ) 1128 { 1129 old = new; 1130 } else 1131 { 1132 if ( kind != 0 ) 1133 { 1134 int io; 1135 for ( int in = 1; in <= size( new[1] ); in ++ ) 1136 { 1137 if( size( new[1][in][2] + old[2] ) > size( old[2] ) ) 1138 { 1139 old[1] = old[1] + list(new[1][in]); 1140 old[2] = old[2] + new[1][in][2]; 1141 } else 1142 { 1143 for( io = 1; io <= size( old[1] ); io ++ ) 1144 { 1145 if( old[1][io][2] == new[1][in][2] ) 1146 { 1147 old[1][io][3] = mergePolysUniq( old[1][io][3], new[1][in][3] ); 1148 break; 1149 }; 1150 }; 1151 }; 1152 }; 1153 } else 1154 { 1155 old[1] = old[1] + new[1]; 1156 old[2] = old[2] + new[2]; 1157 }; 1158 }; 1159 }; 1160 1161 ECall ("mergeRecordsUniq", "{new leads}"); 1162 return (old); 1163 }; 1164 1165 1166 1167 ////////////////////////////////////////////////////////////////////////// 1082 }; 1083 }; 1084 1085 ECall ("mergeRecordsUniq", "{new leads}"); 1086 return (old); 1087 }; 1088 1089 1090 1091 /******************************************************/ 1168 1092 static proc init_bads(int @deg) 1169 1093 " … … 1171 1095 " 1172 1096 { 1173 1174 1175 1176 1177 /// new items:1178 /// {1179 /// list of leads,  empty list1180 /// sum poly "index",  poly(0)1181 /// }1182 1183 1184 }; 1185 1186 / /////////////////////////////////////////////////////////////////////////1097 list @l = list(); 1098 for( int @i = 0; @i <= @deg; @i ++ ) 1099 { 1100 @l[index(@i)] = list( list() , 0); 1101 // new items: 1102 // { 1103 // list of leads,  empty list 1104 // sum poly "index",  poly(0) 1105 // } 1106 }; 1107 return (@l); 1108 }; 1109 1110 /******************************************************/ 1187 1111 static proc zAddBad ( list @leads, list @newz, int @maxDeg, int kind ) 1188 "1112 " 1189 1113 input: 1190 1114 @leads: graded by deg list of … … 1211 1135 " 1212 1136 { 1213 BCall ("zAddBad", "{old leads}, { " + string(@newz) + " }, " + string(@maxDeg) + ", " + string(kind) + ""); 1214 1215 int @newSize = size(@newz); 1216 if ( @newSize < 1 ) 1217 { 1218 return (@leads); 1219 }; 1220 1221 int @i, @j, @k, @dd, @deg; 1222 1223 1224 poly @m, @tm, @ttm; 1225 intvec @exp, @texp, @ttexp; 1226 1227 list @newzz; 1228 int @size, @newdeg, @mydeg; 1229 1230 poly @sum_old, @sum_new; 1231 1232 poly a, b, @temp; 1233 1234 /* 1235 if kind = 0 => for PBW  no products 1236 if kind = 1 => for zReduce  with products 1237 */ 1238 1239 poly pr = 0; int pri; list prs; 1240 1241 for ( @i = @newSize; @i > 0; @i  ) 1242 {// for every new poly (@newz[@i]) do 1243 1244 @m = leadmonom( @newz[@i] ); 1245 @deg = deg(@m); 1246 @exp = leadexp(@m); 1247 1248 //// newzz void new list 1249 @newzz = init_bads(@maxDeg ); 1250 1251 for( @mydeg = @deg; @mydeg <= @maxDeg; @mydeg = @mydeg + @deg ) 1252 {/// adding all possiblities for @newz[@i]^@j; 1253 1254 if ( @mydeg > @deg ) 1137 BCall ("zAddBad", "{old leads}, { " + string(@newz) + " }, " + string(@maxDeg) + ", " + string(kind) + ""); 1138 1139 int @newSize = size(@newz); 1140 if ( @newSize < 1 ) 1141 { 1142 return (@leads); 1143 }; 1144 1145 int @i, @j, @k, @dd, @deg; 1146 1147 1148 poly @m, @tm, @ttm; 1149 intvec @exp, @texp, @ttexp; 1150 1151 list @newzz; 1152 int @size, @newdeg, @mydeg; 1153 1154 poly @sum_old, @sum_new; 1155 1156 poly a, b, @temp; 1157 1158 /* 1159 if kind = 0 => for PBW  no products 1160 if kind = 1 => for zReduce  with products 1161 */ 1162 1163 poly pr = 0; int pri; list prs; 1164 1165 for ( @i = @newSize; @i > 0; @i  ) 1166 {// for every new poly (@newz[@i]) do 1167 1168 @m = leadmonom( @newz[@i] ); 1169 @deg = deg(@m); 1170 @exp = leadexp(@m); 1171 1172 // newzz void new list 1173 @newzz = init_bads(@maxDeg ); 1174 1175 for( @mydeg = @deg; @mydeg <= @maxDeg; @mydeg = @mydeg + @deg ) 1176 {// adding all possiblities for @newz[@i]^@j; 1177 1178 if ( @mydeg > @deg ) 1179 { 1180 @texp = @texp + @exp; 1181 @tm = LM ( @texp ); 1182 if ( kind != 0) 1183 { 1184 pr = QNF_poly( pr * @newz[@i] ); // degrees must be there!!! 1185 }; 1186 } else 1187 { 1188 @texp = @exp; 1189 @tm = @m; 1190 if ( kind != 0) 1191 { 1192 pr = @newz[@i]; 1193 }; 1194 }; 1195 1196 @temp = QNF_poly(@tm); 1197 if( @temp != @tm ) 1198 { 1199 break; 1200 }; 1201 1202 /*!!*/ @newzz[index(@mydeg)] = 1203 /*!!*/ addRecordUniq( @leads[index(@mydeg)], @newzz[index(@mydeg)], @texp, @tm, kind, list(pr) ); 1204 1205 for ( @dd = 1; (@dd <= @maxDeg) and ((@dd + @mydeg) <= @maxDeg ); @dd ++ ) 1206 { // for every good "deg" 1207 1208 @newdeg = @dd + @mydeg; // any deg should be additive!!! 1209 1210 for ( @k = size(@leads[index(@dd)][1]); @k > 0; @k  ) 1211 { 1212 1213 @ttexp = (@leads[index(@dd)][1][@k][1]) + @texp; 1214 @ttm = LM (@ttexp); 1215 1216 if ( kind != 0 ) 1217 { 1218 prs = list(); 1219 1220 for( pri = 1; pri <= size(@leads[index(@dd)][1][@k][3]); pri++) 1255 1221 { 1256 @texp = @texp + @exp; 1257 @tm = LM ( @texp ); 1258 if ( kind != 0) 1259 { 1260 pr = QNF_poly( pr * @newz[@i] ); /// degrees must be there!!! 1261 }; 1262 } else 1263 { 1264 @texp = @exp; 1265 @tm = @m; 1266 if ( kind != 0) 1267 { 1268 pr = @newz[@i]; 1269 }; 1222 // to do products into list and add one list !!! 1223 a = QNF_poly( pr*@leads[index(@dd)][1][@k][3][pri]); 1224 b = QNF_poly( @leads[index(@dd)][1][@k][3][pri]*pr); 1225 1226 prs= prs + list(a); 1227 1228 if ( a!=b ) 1229 { 1230 prs= prs + list(b); 1231 }; 1270 1232 }; 1271 1233 1272 @temp = QNF_poly(@tm); 1273 if( @temp != @tm ) 1274 { 1275 break; 1276 }; 1277 1278 /*!!*/ @newzz[index(@mydeg)] = 1279 /*!!*/ addRecordUniq( @leads[index(@mydeg)], @newzz[index(@mydeg)], @texp, @tm, kind, list(pr) ); 1280 1281 for ( @dd = 1; (@dd <= @maxDeg) and ((@dd + @mydeg) <= @maxDeg ); @dd ++ ) 1282 { //// for every good "deg" 1283 1284 @newdeg = @dd + @mydeg; /// any deg should be additive!!! 1285 1286 for ( @k = size(@leads[index(@dd)][1]); @k > 0; @k  ) 1287 { 1288 1289 @ttexp = (@leads[index(@dd)][1][@k][1]) + @texp; 1290 @ttm = LM (@ttexp); 1291 1292 if ( kind != 0 ) 1293 { 1294 prs = list(); 1295 1296 for( pri = 1; pri <= size(@leads[index(@dd)][1][@k][3]); pri++) 1297 { 1298 /// to do products into list and add one list !!! 1299 a = QNF_poly( pr*@leads[index(@dd)][1][@k][3][pri]); 1300 b = QNF_poly( @leads[index(@dd)][1][@k][3][pri]*pr); 1301 1302 prs= prs + list(a); 1303 1304 if ( a!=b ) 1305 { 1306 prs= prs + list(b); 1307 }; 1308 }; 1309 1310 } else 1311 { 1312 prs = list(pr); 1313 }; 1314 1315 @temp = QNF_poly(@ttm); 1316 if( @temp == @ttm ) 1317 { 1318 1319 /*!!*/ @newzz[index(@newdeg)] = 1320 /*!!*/ addRecordUniq( @leads[index(@newdeg)], @newzz[index(@newdeg)], @ttexp, @ttm, kind, prs ); 1321 }; 1322 1323 1324 }; 1325 }; /// for 1326 1327 if ( @deg == 0 ) 1328 { 1329 break; 1330 }; 1331 }; 1332 1333 for ( @mydeg = 0; @mydeg <= @maxDeg; @mydeg ++ ) 1334 { /// adding currently generated to result 1335 @leads[index(@mydeg)] = mergeRecordsUniq ( @leads[index(@mydeg)], @newzz[index(@mydeg)], kind ); 1336 }; 1337 1338 }; 1339 1340 ECall ("zAddBad", "{new leads}"); 1341 1342 return (@leads); 1343 }; 1344 1345 1346 ////////////////////////////////////////////////////////////////////////// 1347 /// procedure for reducing a given poly wrt already calculated "badleadings" 1348 ////////////////////////////////////////////////////////////////////////// 1349 1350 1351 ////////////////////////////////////////////////////////////////////////// 1234 } else 1235 { 1236 prs = list(pr); 1237 }; 1238 1239 @temp = QNF_poly(@ttm); 1240 if( @temp == @ttm ) 1241 { 1242 1243 /*!!*/ @newzz[index(@newdeg)] = 1244 /*!!*/ addRecordUniq( @leads[index(@newdeg)], @newzz[index(@newdeg)], @ttexp, @ttm, kind, prs ); 1245 }; 1246 1247 1248 }; 1249 }; // for 1250 1251 if ( @deg == 0 ) 1252 { 1253 break; 1254 }; 1255 }; 1256 1257 for ( @mydeg = 0; @mydeg <= @maxDeg; @mydeg ++ ) 1258 { // adding currently generated to result 1259 @leads[index(@mydeg)] = mergeRecordsUniq ( @leads[index(@mydeg)], @newzz[index(@mydeg)], kind ); 1260 }; 1261 1262 }; 1263 1264 ECall ("zAddBad", "{new leads}"); 1265 1266 return (@leads); 1267 }; 1268 1269 1270 /******************************************************/ 1271 // procedure for reducing a given poly wrt already calculated "badleadings" 1272 1273 /******************************************************/ 1352 1274 static proc zReducePoly ( list leads, poly new, poly anfang ) 1353 "1275 " 1354 1276 reduce poly new wrt found leads, 1355 1277 return: list of all possible reductions... 1356 1278 " 1357 1279 { 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1280 poly temp = new; 1281 poly rest = anfang; 1282 1283 poly lm; 1284 number lc; 1285 1286 list LEADS; 1287 1288 int i, n; 1289 list prs; 1290 1291 while ( size(temp) > 0 ) 1292 { 1293 // do for every term... beginning with leading... 'till smallest 1294 lm = leadmonom( temp ); 1295 1296 LEADS = leads[index(deg( lm ))]; // currently users bad leading products... 1297 1298 if( size(LEADS[2] + lm ) <= size(LEADS[2]) ) 1299 { // need to reduce, since leacmonom already in LEADS 1378 1300 1379 for ( i = 1; i <= size(LEADS[1]); i++ ) 1380 { 1381 if( LEADS[1][i][2] == lm ) 1301 for ( i = 1; i <= size(LEADS[1]); i++ ) 1302 { 1303 if( LEADS[1][i][2] == lm ) 1304 { 1305 1306 lc = leadcoef( temp ); // no need be the unit 1307 1308 prs = LEADS[1][i][3]; // shouldbe generated by zAddBad with kind == 1 1309 n = size(prs) ; 1310 1311 if ( n == 1 ) 1312 { // no recursion 1313 1314 temp = leadcoef(prs[1]) * temp  lc * prs[1]; // leadmonom goes down 1315 1316 } else 1317 { // do recursion 1318 1319 list result = list(); 1320 poly newnew; 1321 int f_addrest = 0; 1322 1323 for( int pri = 1; pri <= n ; pri ++ ) 1324 { 1325 newnew = leadcoef(prs[pri]) * temp  lc * prs[pri]; // leadmonom goes down 1326 1327 if( size( newnew ) > 0 ) 1328 { 1329 result = result + zReducePoly(leads, newnew, rest); 1330 } else 1382 1331 { 1383 1384 lc = leadcoef( temp ); // no need be the unit 1385 1386 prs = LEADS[1][i][3]; // shouldbe generated by zAddBad with kind == 1 1387 n = size(prs) ; 1388 1389 if ( n == 1 ) 1390 { // no recursion 1391 1392 temp = leadcoef(prs[1]) * temp  lc * prs[1]; // leadmonom goes down 1393 1394 } else 1395 { // do recursion 1396 1397 list result = list(); 1398 poly newnew; 1399 int f_addrest = 0; 1400 1401 for( int pri = 1; pri <= n ; pri ++ ) 1402 { 1403 newnew = leadcoef(prs[pri]) * temp  lc * prs[pri]; // leadmonom goes down 1404 1405 if( size( newnew ) > 0 ) 1406 { 1407 result = result + zReducePoly(leads, newnew, rest); 1408 } else 1409 { 1410 f_addrest = 1; 1411 }; 1412 }; 1413 1414 if ( f_addrest == 1 ) 1415 { 1416 result = result + list(rest); 1417 }; 1418 return ( result ); 1419 1420 }; 1421 break; 1332 f_addrest = 1; 1422 1333 }; 1423 }; 1424 1425 } else 1426 { // no such leadmonom in leads 1427 1428 rest = rest + lead ( temp ); 1429 temp = temp  lead ( temp ); // leadcoeff goes down 1334 }; 1335 1336 if ( f_addrest == 1 ) 1337 { 1338 result = result + list(rest); 1339 }; 1340 return ( result ); 1341 1342 }; 1343 break; 1430 1344 }; 1431 }; 1432 1433 return (list(rest)); 1434 1435 }; 1436 1437 1438 ////////////////////////////////////////////////////////////////////////// 1345 }; 1346 1347 } else 1348 { // no such leadmonom in leads 1349 1350 rest = rest + lead ( temp ); 1351 temp = temp  lead ( temp ); // leadcoeff goes down 1352 }; 1353 }; 1354 1355 return (list(rest)); 1356 1357 }; 1358 1359 1360 /******************************************************/ 1439 1361 static proc zCancel ( list @tt, list @leads ) 1440 "1362 " 1441 1363 just kill entries of plane PBW base with leading monomials in @leads... 1442 1364 " 1443 1365 { 1444 int @i; 1445 1446 if ( (size(@tt) == 0) ) 1447 { 1448 return (@tt); 1449 }; 1450 1451 list result = list(); 1452 poly g, f; 1453 1454 for ( @i = size(@tt); @i > 0 ; @i  ) 1455 /// for all PBW entries: 1456 { 1457 g = leadmonom(@tt[@i][1]); /// pbw entry 1458 f = @leads[index(deg(g))][2]; /// 'index' poly from @leads 1459 1460 if ( size(f + g) > size(f) ) /// if g not in @leads (works only for monomials) 1461 { 1462 result = list( @tt[@i] ) + result; 1463 }; 1464 }; 1465 1466 return (result); 1467 1468 }; 1469 1470 ////////////////////////////////////////////////////////////////////////// 1471 /// procedure for computing "minimal" set of generators... 1472 ////////////////////////////////////////////////////////////////////////// 1473 1474 ////////////////////////////////////////////////////////////////////////// 1366 int @i; 1367 1368 if ( (size(@tt) == 0) ) 1369 { 1370 return (@tt); 1371 }; 1372 1373 list result = list(); 1374 poly g, f; 1375 1376 for ( @i = size(@tt); @i > 0 ; @i  ) 1377 // for all PBW entries: 1378 { 1379 g = leadmonom(@tt[@i][1]); // pbw entry 1380 f = @leads[index(deg(g))][2]; // 'index' poly from @leads 1381 1382 if ( size(f + g) > size(f) ) // if g not in @leads (works only for monomials) 1383 { 1384 result = list( @tt[@i] ) + result; 1385 }; 1386 }; 1387 1388 return (result); 1389 1390 }; 1391 1392 /******************************************************/ 1393 // procedure for computing "minimal" set of generators... 1394 1395 /******************************************************/ 1475 1396 static proc zReduce( list @z ) 1476 "1397 " 1477 1398 reduce a set @z  base of Center as V.S. 1478 1399 into a minimal set!! 1479 1400 " 1480 1401 { 1481 BCall( "zReduce", @z ); 1482 1483 @z = zRefine ( @z ); 1484 int n = size( @z ); 1485 1486 if ( n < 2 ) 1487 { 1488 return (@z); 1489 }; 1490 1491 int d = maxDeg( stdDeg, @z ); 1492 1493 list leads = init_bads( d ); 1494 1495 list @red; 1496 1497 int f_add; int j; 1498 1499 list result = list(); 1500 1501 poly p; 1502 1503 for ( int i = 1; i <= n ; i++ ) 1504 {// in this order... from least to maximal... 1505 1506 p = @z[i]; 1507 1508 @red = zReducePoly ( leads, p, 0 ); 1509 1510 f_add = 1; 1402 BCall( "zReduce", @z ); 1403 1404 @z = zRefine ( @z ); 1405 int n = size( @z ); 1406 1407 if ( n < 2 ) 1408 { 1409 return (@z); 1410 }; 1411 1412 int d = maxDeg( @z ); 1413 1414 if( d== 1 ) 1415 { 1416 return (@z); 1417 }; 1418 1419 list leads = init_bads( d ); 1420 1421 list @red; 1422 1423 int f_add; int j; 1424 1425 list result = list(); 1426 1427 poly p; 1428 1429 for ( int i = 1; i <= n ; i++ ) 1430 {// in this order... from least to maximal... 1431 1432 p = @z[i]; 1433 1434 @red = zReducePoly ( leads, p, 0 ); 1435 1436 f_add = 1; 1511 1437 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 @red = zRefine( @red ); /// there will be no zeroes after ???1522 1523 1524 1525 1526 1527 1528 1529 {// we got something new....1530 result = result + list(@red); // ??? which one to add? => trying to add all1531 };1532 }; 1533 1534 ECall( "zReduce", result );1535 1536 return (result);1537 };1538 1539 1540 //////////////////////////////////////////////////////////////////////////// 1541 / ///////////////////////////////////////////////////////////////////////////1438 for( j = 1; j <= size(@red); j++ ) 1439 { 1440 if ( @red[j] == 0 ) 1441 { 1442 f_add = 0; 1443 break; // nothing new.... 1444 }; 1445 }; 1446 1447 @red = zRefine( @red ); // there will be no zeroes after ??? 1448 1449 if( size(@red) > 0 ) 1450 { 1451 leads = zAddBad( leads, @red, d, 1); 1452 }; 1453 1454 if ( f_add == 1 ) 1455 { 1456 // we got something new.... 1457 result = result + list(@red); // ??? which one to add? => trying to add all 1458 }; 1459 }; 1460 1461 ECall( "zReduce", result ); 1462 1463 return (result); 1464 }; 1465 1466 1467 /******************************************************/ 1542 1468 // PBW procedures .... 1543 1469 … … 1547 1473 */ 1548 1474 1549 1550 ///////////////////////////////////////////////////////////////////////////// 1475 /******************************************************/ 1551 1476 static proc PBW_init() 1552 "1477 " 1553 1478 PBW[ index(0) ] = PBW part of degree 0: 1554 1479 record: … … 1558 1483 " 1559 1484 { 1560 1561 1562 }; 1563 1564 1565 / ////////////////////////////////////////////////////////////////////////////1485 BCall( "PBW_init" ); 1486 return (list(list(1, 0, 1))); 1487 }; 1488 1489 1490 /******************************************************/ 1566 1491 static proc PBW_sys(list VARS, poly p) 1567 "1492 " 1568 1493 calculate the array [1..NVARS] of records: 1569 1494 record[i] = … … 1576 1501 { 1577 1502 1578 1579 1580 1581 1582 1583 1584 1585 1503 BCall( "PBW_sys" ); 1504 1505 poly t; 1506 for (int v = size(VARS); v>0; v  ) 1507 { 1508 t = VARS[v]; 1509 VARS[v] = list( t, QNF_poly( p*t  t*p ), deg(t) ) ; 1510 }; 1586 1511 1587 1588 }; 1589 1590 / ////////////////////////////////////////////////////////////////////////////1512 return (VARS); 1513 }; 1514 1515 /******************************************************/ 1591 1516 static proc PBW_done( list PBW ) 1592 "1517 " 1593 1518 collects all together, from graded lists into plane list. 1594 1519 … … 1596 1521 " 1597 1522 { 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 }; 1610 1611 / ////////////////////////////////////////////////////////////////////////////1523 BCall( "PBW_done" ); 1524 1525 list result = list(); 1526 1527 int n = size(PBW); 1528 for ( int i = 1; i <= n; i++) 1529 { 1530 result = result + PBW[i]; 1531 }; 1532 1533 return (result); 1534 }; 1535 1536 /******************************************************/ 1612 1537 static proc PBW_next ( list PBW, int k, list sys) 1613 1538 { 1614 BCall( "PBW_next", k ); 1615 1616 list temp; 1617 int kk, nn, ii; 1618 list result = list(); // PBW[ index(k) ] ought to be empty ?? 1619 int N = size(sys); 1620 1621 for ( int i = 1; i <= N; i ++ ) // for all vars 1622 { 1623 kk = (k  sys[i][3]); // any deg should be additive function wrt multiplication 1624 if ( kk >= 0 ) 1539 BCall( "PBW_next", k ); 1540 1541 list temp; 1542 int kk, nn, ii; 1543 list result = list(); // PBW[ index(k) ] ought to be empty ?? 1544 int N = size(sys); 1545 1546 for ( int i = 1; i <= N; i ++ ) // for all vars 1547 { 1548 kk = (k  sys[i][3]); // any deg should be additive function wrt multiplication 1549 if ( kk >= 0 ) 1550 { 1551 nn = size( PBW[ index(kk) ] ); 1552 for ( ii = 1; ii <= nn; ii ++ ) 1553 { 1554 temp = PBW[ index(kk) ][ii]; 1555 1556 if ( temp[3] <= i ) 1625 1557 { 1626 nn = size( PBW[ index(kk) ] ); 1627 for ( ii = 1; ii <= nn; ii ++ ) 1628 { 1629 temp = PBW[ index(kk) ][ii]; 1630 1631 if ( temp[3] <= i ) 1632 { 1633 temp[2] = temp[2]*sys[i][1] + temp[1]*sys[i][2]; // recursive [,] 1634 temp[1] = temp[1]*sys[i][1]; 1635 temp[3] = i; 1558 temp[2] = temp[2]*sys[i][1] + temp[1]*sys[i][2]; // recursive [,] 1559 temp[1] = temp[1]*sys[i][1]; 1560 temp[3] = i; 1636 1561 1637 result = result + list ( temp ); 1638 }; 1639 }; 1562 result = result + list ( temp ); 1640 1563 }; 1641 }; 1642 1643 result = QNF_list( result, 2 + 4); 1644 ECall( "PBW_next", "list result" ); 1645 return (result); 1646 1647 }; 1648 1649 ///////////////////////////////////////////////////////////////////////////// 1564 }; 1565 }; 1566 }; 1567 1568 result = QNF_list( result, 2 + 4); 1569 ECall( "PBW_next", "list result" ); 1570 return (result); 1571 1572 }; 1573 1574 /******************************************************/ 1650 1575 static proc PBW_base( int MaxDeg, poly p ) 1651 1576 { 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1577 BCall( "PBW_base", MaxDeg, p ); 1578 1579 list PBW = list(); 1580 list SYS = PBW_sys( my_vars(), p ); 1581 1582 PBW[index(0)] = PBW_init(); 1583 1584 for (int k = 1; k <= MaxDeg; k ++ ) 1585 { 1586 PBW[index(k)] = PBW_next( PBW, k, SYS ); 1587 }; 1588 1589 return (PBW_done( PBW )); 1665 1590 }; 1666 1591 1667 1592 1668 ///////////////////////////////////////////////////////////////////////////// 1669 ///////////////////////////////////////////////////////////////////////////// 1670 // standard center procedures... "zCommon" 1671 1672 1673 /******************************************************************************/ 1674 static proc my_calc (list @given, poly @BASE_POLY, list #) 1593 /******************************************************/ 1594 // standard center procedures... 1595 1596 /******************************************************/ 1597 static proc FSS (list @given, poly @BASE_POLY) 1675 1598 " 1676 1599 Gauss with computation of kernel v.s basis 1677 " 1678 { 1679 list @nums = list (); 1680 intvec @ones; 1681 int @j, @k, 1682 @n, @v, 1683 @a, @nn; 1600 to be optimizes 1601 " 1602 { 1603 BCall( "FSS", "list", "basepoly: ", @BASE_POLY); 1604 1605 list @nums = list (); 1606 intvec @ones; 1607 int @j, @k, 1608 @n, @v, 1609 @a, @nn; 1684 1610 1685 1686 1611 @n = size( @BASE_POLY ); 1612 @v = size( @given ); 1687 1613 1688 if ( @v == 0)1614 if ( (@v == 0) or (@n == 0) ) 1689 1615 { 1690 return ( list() ); 1691 }; 1692 1693 if ( size(#)>0 ) 1694 { 1695 poly @Next_Ad_var; 1696 if ( typeof(#[1]) == "int") 1697 { 1698 int @z_next = #[1]; 1699 @Next_Ad_var = var (@z_next); 1700 }; 1701 1702 if ( typeof(#[1]) == "poly") 1703 { 1704 @Next_Ad_var = #[1]; 1705 }; 1706 }; 1616 return (@given); 1617 }; 1618 1619 matrix MD[@n][@v]; 1620 1621 number @max; 1622 poly @test; 1623 poly @t; 1624 1625 list LM = list(); // rec: { exp, poly } 1626 1627 for( @k = @v; @k > 0; @k  ) 1628 { 1629 LM[@k] = list(); 1630 @t = @given[@k][2]; 1631 LM[@k][2]= @t; 1632 LM[@k][1]= leadexp( @t ); 1633 }; 1634 1635 1636 intvec @w; 1637 for ( @a = 1 ; @a <= @n ; @a ++ ) 1638 { 1639 @w = leadexp(@BASE_POLY[@a]); 1640 for( @k = @v; @k > 0; @k  ) 1641 { 1642 if( @w == LM[@k][1] ) 1643 { 1644 @t = LM[@k][2]; 1645 MD[@a, @k] = leadcoef( @t ); 1646 @t = @t  lead ( @t ); 1647 LM[@k][2] = @t; 1648 1649 LM[@k][1] = leadexp(@t); 1650 }; 1651 }; 1652 1653 }; 1654 1655 int @x, @y; 1656 number @min; 1657 1658 number @div; 1659 // @nums = list(); 1707 1660 1708 poly @t; 1661 // Gauss 1662 // print("Before Gauss, Matrix: "); 1663 // print( MD ); 1664 1709 1665 1710 if ( @n == 0 ) 1711 { 1712 if (defined(@Next_Ad_var)) 1713 { 1714 for ( @k = @v; @k > 0; @k ) 1715 { 1716 @t = @given[@k][1]; 1717 @given[@k][2] = @Next_Ad_var * @t  @t * @Next_Ad_var; 1718 }; 1719 }; 1720 return (@given); 1721 }; 1722 1723 1724 matrix MD[@n][@v]; 1725 1726 number @max; 1727 poly @test; 1728 1729 list LM = list(); /// rec: { exp, poly } 1730 1731 1732 1733 1734 for( @k = @v; @k > 0; @k  ) 1735 { 1736 LM[@k] = list(); 1737 @t = @given[@k][2]; 1738 LM[@k][2]= @t; 1739 LM[@k][1]= leadexp( @t ); 1740 }; 1741 1742 1743 intvec @w; 1744 for ( @a = 1 ; @a <= @n ; @a ++ ) 1745 { 1746 @w = leadexp(@BASE_POLY[@a]); 1747 for( @k = @v; @k > 0; @k  ) 1748 { 1749 if( @w == LM[@k][1] ) 1750 { 1751 @t = LM[@k][2]; 1752 MD[@a, @k] = leadcoef( @t ); 1753 @t = @t  lead ( @t ); 1754 LM[@k][2] = @t; 1755 1756 LM[@k][1] = leadexp(@t); 1757 }; 1758 }; 1759 1760 }; 1761 1762 int @x, @y; 1763 number @min; 1764 1765 number @div; 1766 /// @nums = list(); 1767 1768 /// Gauss 1769 1770 @x = 1; 1771 @y = 1; 1772 while ( (@y <= @n) && (@x <= @v)) 1666 @x = 1; 1667 @y = 1; 1668 while ( (@y <= @n) && (@x <= @v)) 1773 1669 // main Gauss loop... 1774 1670 { 1775 @min = leadcoef( MD[@y, @x] ); /// curr diag.1776 1777 1671 @min = leadcoef( MD[@y, @x] ); // curr diag. 1672 1673 if ( @min == 0 ) 1778 1674 // if zero on diag... 1779 1675 { 1780 1781 1782 /// let's find the minimal1783 1784 { 1785 1786 1676 @j = 0; 1677 1678 // let's find the minimal 1679 for ( @k = @y+1; @k <= @n; @k ++ ) 1680 { 1681 @max = leadcoef( MD[ @k, @x ] ); 1682 if ( @max != 0 ) 1787 1683 { 1788 @j = @k; 1789 @min = @max; 1790 ////// @k ++; 1791 break; //// this pure for 1792 /// continue; // for find min 1684 @j = @k; 1685 @min = @max; 1686 // @k ++; 1687 break; // this pure for 1688 // continue; 1689 // for find min 1793 1690 }; 1794 }; // /for (@k) // found minimal1795 1796 1691 }; // for (@k) // found minimal 1692 1693 if ( @j == 0 ) 1797 1694 { 1798 /// das ist gut! ///1799 1800 1801 1695 // das ist gut! // 1696 @nums = @nums + list ( list ( @x, @y ) ); 1697 @x ++; 1698 continue; // while 1802 1699 } ; 1803 1700 1804 for ( @k = @x ; @k <= @v ; @k ++ ) 1805 { 1806 @t = MD[ @j, @k ]; 1807 MD[ @j, @k ] = MD[ @y, @k ]; 1808 MD[ @y, @k ] = @t; 1809 }; 1701 for ( @k = @x ; @k <= @v ; @k ++ ) 1702 { 1703 @t = MD[ @j, @k ]; 1704 MD[ @j, @k ] = MD[ @y, @k ]; 1705 MD[ @y, @k ] = @t; 1706 }; 1707 1708 }; // if diag === zero. 1709 @ones[@y] = @x; 1810 1710 1811 }; /// if diag === zero. 1812 @ones[@y] = @x; 1813 1814 if ( @min == 0 ) 1815 { 1816 // write_latex_str ( " ******************* ERROR ******************* " ); 1817 quit; 1818 }; 1819 1711 if ( @min == 0 ) 1712 { 1713 // write_latex_str ( " ******************* ERROR ******************* " ); 1714 quit; 1715 }; 1716 1717 1718 if ( @min != 1) // let's norm the row... to make the '1' ! 1719 { 1720 @min = 1 / @min; 1721 for ( @k = @x; @k <= @v; @k++ ) 1722 { 1723 @max = leadcoef( MD[@y, @k] ); 1724 1725 if ( @max == 0) 1726 { 1727 @k ++ ; 1728 continue; // for : norming the row... 1729 }; 1730 1731 MD[@y, @k] = @max * @min; // here must be Field ... 1732 }; 1733 1734 // @min = 1; 1735 1736 }; 1737 1738 // here must be @min != 0; 1739 for ( @k = 1; @k <= @n; @k++ ) 1740 { 1741 @max = leadcoef( MD[@k, @x] ); 1742 1743 if ( ( @k == @y)  (@max == 0) ) 1744 { 1745 @k ++ ; 1746 continue; 1747 }; 1748 1749 for ( @j = @x; @j <= @v ; @j ++ ) 1750 { 1751 MD[ @k, @j ] = MD[ @k, @j ]  @max * MD[ @y, @j ]; 1752 }; 1753 1754 }; //killing 1755 1756 @x ++; 1757 @y ++; 1758 }; //main while. 1759 /******************************************************/ 1820 1760 1821 if ( @min != 1) /// let's norm the row... to make the '1' ! 1822 { 1823 @min = 1 / @min; 1824 for ( @k = @x; @k <= @v; @k++ ) 1825 { 1826 @max = leadcoef( MD[@y, @k] ); 1827 1828 if ( @max == 0) 1829 { 1830 @k ++ ; 1831 continue; /// for : norming the row... 1832 }; 1833 1834 MD[@y, @k] = @max * @min; ////here must be Field ... 1835 }; 1836 1837 //// @min = 1; 1838 1839 }; 1840 1841 //// here must be @min != 0; 1842 for ( @k = 1; @k <= @n; @k++ ) 1843 { 1844 @max = leadcoef( MD[@k, @x] ); 1845 1846 if ( ( @k == @y)  (@max == 0) ) 1847 { 1848 @k ++ ; 1849 continue; 1850 }; 1851 1852 for ( @j = @x; @j <= @v ; @j ++ ) 1853 { 1854 MD[ @k, @j ] = MD[ @k, @j ]  @max * MD[ @y, @j ]; 1855 }; 1856 1857 }; ///killing 1858 1859 @x ++; 1860 @y ++; 1861 }; //////main while. 1862 ////////////////////////////////////////////////////// 1863 1864 1865 /// computation of kernel's basis 1866 1867 1868 if ( @x <= @v ) 1869 { 1870 for ( @k = @x; @k <= @v ; @k ++ ) 1871 { 1872 @nums = @nums + list ( list ( @k, @n+1 ) ); 1761 // print("Gaussed Matrix: "); 1762 // print( MD ); 1763 1764 // computation of kernel's basis 1765 1766 if ( @x <= @v ) 1767 { 1768 for ( @k = @x; @k <= @v ; @k ++ ) 1769 { 1770 @nums = @nums + list ( list ( @k, @n+1 ) ); 1873 1771 } 1874 1772 } 1875 1773 1876 if ( @y <= @n ) 1877 { 1878 @n = @y  1; 1879 }; 1774 // print("Nums: " ); 1775 // print (@nums); 1880 1776 1881 1882 list result = list(); 1777 list result = list(); 1883 1778 1884 /// real calculations of the Base of a Ker as V.S. 1885 1886 for ( @k = 1; @k <= size(@nums) ; @k ++ ) 1887 { 1888 @x = @nums[@k][1]; 1889 @j = @nums[@k][2]; 1890 1891 @t = @given[@x][1]; 1892 1893 for ( @y = 1; @y < @j ; @y ++ ) 1894 /// for every "@x" column 1895 { 1896 @max = leadcoef( MD[@y, @x] ); 1897 if ( (@max != 0) ) 1898 { 1899 @a = @ones[@y]; 1900 @t = @t  @max * @given[@a][1]; 1901 }; 1902 }; 1903 1904 if ( defined(@Next_Ad_var) ) 1905 { 1906 result[@k] = list ( @t, @Next_Ad_var * @t  @t * @Next_Ad_var ) ; 1779 // real calculations of the Base of a Ker as V.S. 1780 1781 for ( @k = 1; @k <= size(@nums) ; @k ++ ) 1782 { 1783 @x = @nums[@k][1]; 1784 @j = @nums[@k][2]; 1785 1786 @t = @given[@x][1]; 1787 1788 for ( @y = 1; @y < @j ; @y ++ ) 1789 // for every "@x" column 1790 { 1791 @max = leadcoef( MD[@y, @x] ); 1792 if ( (@max != 0) ) 1793 { 1794 @a = @ones[@y]; 1795 @t = @t  @max * @given[@a][1]; 1796 }; 1797 }; 1798 1799 result[@k] = @t; 1800 }; 1801 1802 ECall( "FSS", "list" ); 1803 return ( result ); 1804 }; 1805 1806 /******************************************************/ 1807 static proc reduce_one_per_row(list @given) 1808 { 1809 BCall( "reduce_one_per_row" ); 1810 1811 int @k; 1812 int @l = size (@given); 1813 1814 if( @l == 0 ) 1815 { 1816 return (@given); 1817 }; 1818 1819 int @char = char(basering); 1820 1821 intvec @marks; 1822 1823 if ( @char == 0 ) 1824 { 1825 poly @t = poly(0); 1826 }; 1827 1828 list @unis = list (); 1829 poly p; 1830 1831 for ( @k = @l; @k > 0; @k  ) 1832 { 1833 p = uni_poly( @given[@k][2] ); 1834 @unis[@k] = p; 1835 1836 if( p == 0 ) 1837 { 1838 @marks[@k] = 2; 1907 1839 } else 1908 { 1909 result[@k] = @t; 1910 } 1911 }; 1912 1913 return ( result ); 1914 }; 1915 1916 1917 /******************************************************************************/ 1918 /* 1919 static proc add2list( poly l, poly q ) 1920 " 1921 l is a sorted list of monomials (polynom with coefs: 0 or 1) 1922 adding monomials from q, checking also for presence 1923 1924 ?? NEED DEBUG ?? 1925 " 1926 { 1927 BCall( "add2list", "{" + string(l) + "} , " + string(q) ); 1928 1929 poly t; 1930 for( int i = size(q); i>0; i ) 1931 { 1932 1933 /// trick: 1934 /// check for presens of a monomial q[i] in l 1935 /// if char = 0 then new_size > (if not present) or = (if already present) 1936 /// if char > 0 then new_size > (if not present) or =< (if already present) 1937 t = l + leadmonom(q[i]); 1938 if ( size(t) > size(l) ) 1840 { 1841 @marks[@k] = 1; 1842 if ( @char == 0 ) 1843 { 1844 @t = @t + p; 1845 }; 1846 }; 1847 }; 1848 1849 if ( @char != 0 ) 1850 { 1851 def save = basering; 1852 execute( "ring NewRingWithGoodField = (0), (" + varstr(save) + "), (" + ordstr(save) + "); "); 1853 poly @t = 0; 1854 1855 if(! defined (@unis) ) 1856 { 1857 list @unis = imap( save, @unis ); 1858 }; 1859 1860 for ( @k = @l; @k > 0; @k  ) 1861 { 1862 if( @marks[@k] == 1 ) 1863 { 1864 @t = @t + @unis[@k]; 1865 }; 1866 }; 1867 }; 1868 1869 int @loop_size = size(@t); 1870 poly @for_delete, @tt; 1871 int @ll; 1872 1873 while( @loop_size > 0 ) 1874 { 1875 @for_delete = poly(0); 1876 @ll = size(@t); 1877 1878 for ( @k = @ll; @k > 0; @k  ) 1879 { 1880 if ( leadcoef(@t[@k]) == 1 ) 1881 { 1882 @for_delete = @for_delete + @t[@k]; 1883 }; 1884 }; 1885 1886 @loop_size = size( @for_delete ); 1887 1888 if ( @loop_size>0 ) 1889 { 1890 for( @k = @l ; @k > 0 ; @k  ) 1891 { 1892 if ( @marks[@k] == 1) 1939 1893 { 1940 l = t; 1941 }; 1942 }; 1943 1944 ECall( "add2list", "{" + string(l) + "}" ); 1945 return (l); 1946 }; 1947 */ 1948 1949 /******************************************************************************/ 1950 /* 1951 static proc rm4list( poly l, poly q ) 1952 " 1953 l is a sorted list of monomials (polynom with coefs: 0 or 1) 1954 removing from l monomials occuring in q 1955 1956 ?? NEED DEBUG ?? 1957 " 1958 { 1959 BCall( "rm4list", "{" + string(l) + "} , " + string(q) ); 1960 poly t; 1961 for( int i = size(q); i>0; i ) 1962 { 1963 1964 /// trick: 1965 /// check for presens of a monomial q[i] in l 1966 /// if char = 0 then new_size > (if not present) or = (if already present) 1967 /// if char > 0 then new_size > (if not present) or =< (if already present) 1968 t = leadmonom(q[i]); 1969 if ( size(l+t) <= size(l) ) 1970 { 1971 l = t  t; 1972 if( size(l) == 0 ) 1973 { 1974 break; 1975 }; 1976 }; 1977 }; 1978 1979 ECall( "rm4list", "{" + string(l) + "}" ); 1980 return (l); 1981 }; 1982 */ 1983 1984 /******************************************************************************/ 1985 static proc reduce_one_per_row(list @given) 1986 { 1987 BCall( "reduce_one_per_row" ); 1988 1989 int @k; 1990 int @l = size (@given); 1991 1992 if( @l == 0 ) 1993 { 1994 return (@given); 1995 }; 1996 1997 int @char = char(basering); 1998 1999 intvec @marks; 2000 2001 if ( @char == 0 ) 2002 { 2003 poly @t = poly(0); 2004 }; 2005 2006 list @unis = list (); 2007 poly p; 2008 2009 for ( @k = @l; @k > 0; @k  ) 2010 { 2011 p = uni_poly( @given[@k][2] ); 2012 @unis[@k] = p; 2013 2014 if( p == 0 ) 2015 { 2016 @marks[@k] = 2; 2017 } else 2018 { 2019 @marks[@k] = 1; 2020 if ( @char == 0 ) 2021 { 2022 @t = @t + p; 2023 }; 1894 @tt = @unis[@k]; 1895 1896 if( size( @for_delete + @tt ) != ( size( @for_delete ) + size( @tt ) ) ) 1897 { 1898 @t = @t  @tt; 1899 @marks[@k] = 0; 1900 }; 2024 1901 }; 2025 }; 2026 2027 if ( @char != 0 ) 2028 { 2029 def save = basering; 2030 execute( "ring NewRingWithGoodField = (0), (" + varstr(save) + "), (" + ordstr(save) + "); "); 2031 poly @t = 0; 2032 2033 if(! defined (@unis) ) 2034 { 2035 list @unis = imap( save, @unis ); 2036 }; 2037 2038 for ( @k = @l; @k > 0; @k  ) 2039 { 2040 if( @marks[@k] == 1 ) 2041 { 2042 @t = @t + @unis[@k]; 2043 }; 2044 }; 2045 }; 2046 2047 int @loop_size = size(@t); 2048 poly @for_delete, @tt; 2049 int @ll; 2050 2051 while( @loop_size > 0 ) 2052 { 2053 @for_delete = poly(0); 2054 @ll = size(@t); 2055 2056 for ( @k = @ll; @k > 0; @k  ) 2057 { 2058 if ( leadcoef(@t[@k]) == 1 ) 2059 { 2060 @for_delete = @for_delete + @t[@k]; 2061 }; 2062 }; 2063 2064 @loop_size = size( @for_delete ); 2065 2066 if ( @loop_size>0 ) 2067 { 2068 for( @k = @l ; @k > 0 ; @k  ) 2069 { 2070 if ( @marks[@k] == 1) 2071 { 2072 @tt = @unis[@k]; 2073 2074 if( size( @for_delete + @tt ) != ( size( @for_delete ) + size( @tt ) ) ) 2075 { 2076 @t = @t  @tt; 2077 @marks[@k] = 0; 2078 }; 2079 }; 2080 }; 2081 }; 2082 }; 2083 2084 if ( @char != 0 ) 2085 { 2086 setring(save); 2087 kill(NewRingWithGoodField); 2088 }; 2089 2090 list @reduced = list(); 2091 2092 for ( @k = @l ; @k>0 ; @k ) 2093 { 2094 if (@marks[@k]==2) 2095 { 2096 @reduced = list ( @given[@k] ) + @reduced ; 2097 } else 2098 { 2099 if (@marks[@k]==1) 2100 { 2101 @reduced = @reduced + list ( @given[@k] ); 2102 }; 2103 }; 2104 }; 1902 }; 1903 }; 1904 }; 1905 1906 if ( @char != 0 ) 1907 { 1908 setring(save); 1909 kill(NewRingWithGoodField); 1910 }; 1911 1912 list @reduced = list(); 1913 1914 for ( @k = @l ; @k>0 ; @k ) 1915 { 1916 if (@marks[@k]==2) 1917 { 1918 @reduced = list ( @given[@k] ) + @reduced ; 1919 } else 1920 { 1921 if (@marks[@k]==1) 1922 { 1923 @reduced = @reduced + list ( @given[@k] ); 1924 }; 1925 }; 1926 }; 2105 1927 2106 ECall( "reduce_one_per_row", "{...l...}" );2107 2108 }; 2109 2110 2111 2112 /****************************************************** ************************/2113 static proc calc_base (list @AD_GIVEN , list #)2114 "1928 ECall( "reduce_one_per_row", "structured list" ); 1929 return (@reduced); 1930 }; 1931 1932 1933 1934 /******************************************************/ 1935 static proc calc_base (list @AD_GIVEN) 1936 " 2115 1937 sort out given 'pbw' into groups due to common monoms in images = ([2])s 2116 1938 " 2117 1939 { 2118 BCall( "calc_base" ); 2119 2120 if ( (size(#)>0) ) 2121 { 2122 int @z_next; 2123 poly @z_nextp; 2124 2125 if ( typeof(#[1]) == "int") 2126 { 2127 @z_next = #[1]; 2128 @z_nextp = var(@z_next); 2129 }; 2130 2131 if ( typeof(#[1]) == "poly") 2132 { 2133 @z_nextp = #[1]; 2134 }; 2135 2136 }; 2137 2138 2139 list @AD_CALC = list(); 2140 int @ll, @k, @j, @loop_size; 2141 2142 poly @t = poly(0); 2143 poly @tt, @for_delete; 2144 2145 int @t_size; 2146 2147 2148 list @GR, @GR_TEMP; 2149 2150 int @number = size(@AD_GIVEN); 2151 2152 while ( @number > 0 ) 2153 { 2154 @tt = @AD_GIVEN[@number][2]; 2155 if ( size (@tt) == 0) 2156 { 2157 @t = @AD_GIVEN[@number][1]; 2158 if (defined(@z_nextp)) 2159 { 2160 @tt = @z_nextp * @t  @t * @z_nextp; /// right? (?) 2161 @AD_CALC = @AD_CALC + list ( list ( @t, @tt ) ); 2162 } else 2163 { 2164 @AD_CALC = @AD_CALC + list ( @t ); 2165 }; 2166 1940 BCall( "calc_base" ); 1941 1942 list @AD_CALC = list(); 1943 int @ll, @k, @j, @loop_size; 1944 1945 poly @t = poly(0); 1946 poly @tt, @for_delete; 1947 1948 int @t_size; 1949 1950 list @GR, @GR_TEMP; 1951 1952 int @number = size(@AD_GIVEN); 1953 1954 while ( @number > 0 ) 1955 { 1956 @tt = @AD_GIVEN[@number][2]; 1957 if ( size (@tt) == 0) 1958 { 1959 @AD_CALC = @AD_CALC + list ( @AD_GIVEN[@number][1] ); 2167 1960 } else 2168 {1961 { 2169 1962 @t = uni_poly( @tt ); 2170 1963 @t_size = size(@t); … … 2176 1969 @loop_size = size(@GR); 2177 1970 if ( @loop_size == 0 ) 2178 {1971 { 2179 1972 @GR = list(@GR_TEMP); 2180 } else 2181 { 2182 for ( @k = @loop_size; @k > 0 ; @k  ) 1973 } else 2183 1974 { 2184 @tt = @GR[@k][1]; 2185 if ( size( @t + @tt ) != ( @t_size + size(@tt) ) ) 2186 // whether @tt and @i intersencts? ( will not work in char == 2 !!!) 1975 for ( @k = @loop_size; @k > 0 ; @k  ) 2187 1976 { 2188 2189 if ( char(basering) == 0 ) 1977 @tt = @GR[@k][1]; 1978 if ( size( @t + @tt ) != ( @t_size + size(@tt) ) ) 1979 // whether @tt and @i intersencts? ( will not work in char == 2 !!!) 2190 1980 { 2191 @GR_TEMP[1] = @GR_TEMP[1] + @tt; 2192 } else 2193 { 2194 @GR_TEMP[1] = uni_poly( @GR_TEMP[1] + @tt ); 1981 1982 if ( char(basering) == 0 ) 1983 { 1984 @GR_TEMP[1] = @GR_TEMP[1] + @tt; 1985 } else 1986 { 1987 @GR_TEMP[1] = uni_poly( @GR_TEMP[1] + @tt ); 1988 }; 1989 1990 @GR_TEMP[2] = @GR_TEMP[2] + @GR[@k][2]; 1991 @GR = delete ( @GR, @k ); 2195 1992 }; 2196 2197 @GR_TEMP[2] = @GR_TEMP[2] + @GR[@k][2];2198 @GR = delete ( @GR, @k );2199 1993 }; 1994 @GR = @GR + list(@GR_TEMP); 2200 1995 }; 2201 @GR = @GR + list(@GR_TEMP); 2202 }; 2203 }; 2204 @number ; 2205 }; /// main while 1996 }; 1997 @number ; 1998 }; // main while 2206 1999 2207 list @res; 2208 2209 for ( @k = size(@GR); @k > 0 ; @k  ) 2210 { 2211 if ( size (@GR[@k][2]) > 1 ) /// ! zeroes in AD_CALC so here must be non zero 2212 { 2213 if(defined(@z_nextp)) 2214 { 2215 @res = my_calc ( @GR[@k][2], uni_poly(@GR[@k][1]), @z_nextp); 2216 } else 2217 { 2218 @res = my_calc ( @GR[@k][2], uni_poly(@GR[@k][1])); 2219 } 2220 if ( size (@res) > 0 ) 2221 { 2222 @AD_CALC = @AD_CALC + @res; 2000 list @res; 2001 2002 for ( @k = size(@GR); @k > 0 ; @k  ) 2003 { 2004 if ( size (@GR[@k][2]) > 1 ) // ! zeroes in AD_CALC so here must be non zero 2005 { 2006 @res = FSS ( @GR[@k][2], uni_poly(@GR[@k][1])); 2007 2008 if ( size (@res) > 0 ) 2009 { 2010 @AD_CALC = @AD_CALC + @res; 2223 2011 }; 2224 2012 }; 2225 2013 }; 2226 2014 2227 ECall( "calc_base" ); 2228 return (@AD_CALC); 2229 2230 }; 2231 2232 /******************************************************************************/ 2015 ECall( "calc_base" ); 2016 return (@AD_CALC); 2017 2018 }; 2019 2020 2021 /******************************************************/ 2022 static proc ApplyAd( list l, list # ) 2023 " 2024 Apply Ad_{#} to every element of a list of polynomils 2025 " 2026 { 2027 BCall( "ApplyAd", l, # ); 2028 2029 if( (size(l) == 0) or (size(#) == 0) ) 2030 { 2031 return (l); 2032 }; 2033 2034 poly f, t; 2035 2036 if ( typeof(#[1]) == "int") 2037 { 2038 int next = #[1]; 2039 f = my_var (next); 2040 } else 2041 { 2042 if ( typeof(#[1]) == "poly") 2043 { 2044 f = #[1]; 2045 } else 2046 { 2047 print("Error: cannot differentiate with '" + string(#)+"'"); 2048 return (); 2049 }; 2050 }; 2051 2052 for ( int k = size(l); k > 0; k ) 2053 { 2054 t = l[k]; 2055 l[k] = list( t, f * t  t * f ); 2056 }; 2057 2058 ECall("ApplyAd", l ); 2059 return (l); 2060 }; 2061 2062 2063 /******************************************************/ 2233 2064 static proc calc_k_base (list l, list #) 2234 "2235 calculation of Ker asVector Space2236 " 2237 { 2238 2239 2240 list t = reduce_one_per_row( l, 0); /// optimization (a1)2241 return( QNF_list ( calc_base(t, #), 2) ); /// calculation of groupps (a2) + gauss.2242 2243 }; 2244 2245 /******************************************************************************/ 2246 2247 / //////////////////////////////////////////////////////////////////////////////2065 " 2066 calculation of Kernel of an Ad operator as a KVector Space 2067 " 2068 { 2069 BCall( "calc_k_base", "list, " + string(#) ); 2070 2071 list t = reduce_one_per_row( l, 0); // optimization (a1) 2072 return( QNF_list ( ApplyAd( calc_base(t), # ), 2) ); // calculation of groupps (a2) + gauss. 2073 2074 }; 2075 2076 LIB "poly.lib"; // for content in proc makeIdeal 2077 2078 /******************************************************/ 2248 2079 static proc makeIdeal ( list l ) 2249 "2080 " 2250 2081 return: ideal: where the generators are polynomials from list, without 1 and zeroes 2251 2082 " 2252 2083 { 2253 2254 2255 2256 2257 2258 2259 p = l[i][1]; /// just take the 1st polynom...2260 2261 2262 2263 2084 ideal I; poly p; 2085 2086 for( int i = 1; i <= size(l); i++ ) 2087 { 2088 if ( typeof( l[i] ) == "list" ) 2089 { 2090 p = l[i][1]; // just take the 1st polynom... 2091 } else 2092 { 2093 p = l[i]; 2094 }; 2264 2095 2265 if ( (p != 1) and (p != 0) ) 2266 { 2267 I = I, p; 2268 }; 2269 }; 2270 2271 I = simplify ( I, 2 ); /// no zeroes 2272 2273 return(I); 2274 }; 2275 2276 /******************************************************************************/ 2096 p = cleardenom( p* (1/content(p)) ); 2097 if ( (p != 1) and (p != 0) ) 2098 { 2099 I = I, p; 2100 }; 2101 }; 2102 2103 I = simplify ( I, 2 ); // no zeroes 2104 2105 return(I); 2106 }; 2107 2108 /******************************************************/ 2277 2109 static proc inCenter_poly( poly p ) 2278 "2110 " 2279 2111 if p in center => return 1 2280 2112 otherwise return 0 2281 2113 " 2282 2114 { 2283 2284 2285 { 2286 2287 2288 { 2289 2290 2291 2292 2293 2294 }; 2295 }; 2296 2297 2298 { 2299 2300 }; 2301 2302 }; 2303 2304 /****************************************************** ************************/2115 poly t; 2116 for (int k = nvars(basering); k>0; k ) 2117 { 2118 t = var(k); 2119 if ( QNF_poly(t * p  p * t) != 0 ) 2120 { 2121 if( toprint() ) 2122 { 2123 "POLY: ", string (p), " is NOT in center"; 2124 }; 2125 return (0); 2126 }; 2127 }; 2128 2129 if( toprint() ) 2130 { 2131 "POLY: ", string (p), " is in center"; 2132 }; 2133 return (1); 2134 }; 2135 2136 /******************************************************/ 2305 2137 static proc inCenter_list(def l) 2306 2138 { 2307 2308 { 2309 2310 { 2311 2312 { 2313 2139 for ( int @i = 1; @i <= size(l); @i++ ) 2140 { 2141 if ( typeof(l[@i])=="poly" ) 2142 { 2143 if (! inCenter_poly(l[@i]) ) 2144 { 2145 return(0); 2314 2146 }; 2315 2147 2316 2148 } else 2317 {2149 { 2318 2150 if ( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") ) 2319 {2320 2151 { 2152 if (! inCenter_list(l[@i]) ) 2321 2153 { 2322 2154 return(0); 2323 2155 }; 2324 };2325 };2326 }; 2327 2328 }; 2329 2330 /****************************************************** ************************/2156 }; 2157 }; 2158 }; 2159 return(1); 2160 }; 2161 2162 /******************************************************/ 2331 2163 static proc inCentralizer_poly( poly p, poly f ) 2332 "2164 " 2333 2165 if p in centralizer(f) => return 1 2334 2166 otherwise return 0 2335 2167 " 2336 2168 { 2337 2169 if ( QNF_poly(f * p  p * f) != 0 ) 2338 2170 { 2339 2340 2341 2342 2343 2344 }; 2345 2346 2347 { 2348 2349 }; 2350 2351 }; 2352 2353 /****************************************************** ************************/2171 if( toprint() ) 2172 { 2173 "POLY: ", string (p), " is NOT in centralizer(f)"; 2174 }; 2175 return (0); 2176 }; 2177 2178 if( toprint() ) 2179 { 2180 "POLY: ", string (p), " is in centralizer(f)"; 2181 }; 2182 return (1); 2183 }; 2184 2185 /******************************************************/ 2354 2186 static proc inCentralizer_list( def l, poly f ) 2355 2187 { 2356 2188 for ( int @i = 1; @i <= size(l); @i++ ) 2357 {2358 2359 { 2360 2361 2362 2363 2189 { 2190 if ( typeof(l[@i])=="poly" ) 2191 { 2192 if (! inCentralizer_poly(l[@i], f) ) 2193 { 2194 return(0); 2195 }; 2364 2196 2365 2197 } else 2366 {2198 { 2367 2199 if ( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") ) 2368 {2369 2200 { 2201 if (! inCentralizer_list(l[@i], f) ) 2370 2202 { 2371 2203 return(0); 2372 2204 }; 2373 };2374 };2375 };2205 }; 2206 }; 2207 }; 2376 2208 return(1); 2377 2209 }; 2378 2210 2379 2211 2380 / //////////////////////////////////////////////////////////////////////////////2381 / //////////////////////////////////////////////////////////////////////////////2382 // /main algorithms2383 2384 2385 / //////////////////////////////////////////////////////////////////////////////2386 / //////////////////////////////////////////////////////////////////////////////2387 // /center's computation2388 2389 / //////////////////////////////////////////////////////////////////////////////2212 /******************************************************/ 2213 /******************************************************/ 2214 // main algorithms 2215 2216 2217 /******************************************************/ 2218 /******************************************************/ 2219 // center's computation 2220 2221 /******************************************************/ 2390 2222 /* static */ proc center_min_iterative( int MaxDeg, list # ) 2391 "2223 " 2392 2224 computes the 'minimal' set of central elements (of deg <= MaxDeg) in iterative way 2393 2225 Note: based on calc_k_base, zAddBad, zRefine, zCancel, PBW_* 2394 2226 " 2395 2227 { 2396 2397 2398 2399 int m = ( MaxDeg < 0 ); /// no bound on Degree2400 2401 int MinDeg = 6; /// starting guess for MaxDeg2402 int Delta = 4; /// increment of MaxDeg2403 2404 2405 2406 2407 2408 2409 2410 2411 2412 2413 2414 2415 2416 2417 2418 2419 2420 2421 2422 2423 list @z = list (); /// center list2424 list @l = init_bads( MaxDeg ); /// verbotten loeadexps...2425 2426 2427 2428 2429 2430 2431 2432 2433 2434 2435 2436 @q = zRefine (calc_k_base(@q)); /// new center!2437 2438 2439 2440 2441 @z = @z + @q; /// computed central elements2228 BCall("center_min_iterative", MaxDeg, #); 2229 2230 int n = myInt(#); 2231 int m = ( MaxDeg < 0 ); // no bound on Degree 2232 2233 int MinDeg = 6; // starting guess for MaxDeg 2234 int Delta = 4; // increment of MaxDeg 2235 2236 if( m ) 2237 { 2238 // minimal guess 2239 MaxDeg = MinDeg; 2240 }; 2241 2242 list @q; int @i; int N = nvars(basering); 2243 2244 my_var_init(); // setup for my_var(n) 2245 QNF_init(); 2246 2247 list PBW = list(); 2248 list PBW_PLAIN = list(); 2249 2250 PBW[ index(0) ] = PBW_init(); 2251 2252 list SYS = PBW_sys( my_vars(), my_var(1) ); 2253 2254 2255 list @z = list (); // center list 2256 list @l = init_bads( MaxDeg ); // verbotten loeadexps... 2257 2258 @q = PBW[ index(0) ]; 2259 2260 int k = 1; 2261 while( k <= ( MaxDeg+1 ) ) 2262 { 2263 for ( @i = 2; @i <= N; @i ++ ) 2264 { 2265 @q = calc_k_base (@q, my_var(@i)); 2266 }; 2267 2268 @q = zRefine (calc_k_base(@q)); // new center! 2269 2270 2271 if ( size(@q) > 0 ) 2272 { 2273 @z = @z + @q; // computed central elements 2442 2274 2443 2444 2445 2446 2447 2448 2449 2450 2451 2452 2453 break; /// that's all2454 2275 if( (n > 0) and (size(@z) > n) ) 2276 { 2277 break; // found all central elements 2278 }; 2279 }; 2280 2281 if( k == ( MaxDeg+1 ) ) 2282 { 2283 if( (n == 0) and ( !m ) ) 2284 { 2285 break; // that's all 2286 }; 2455 2287 2456 2288 MaxDeg = MaxDeg + Delta; 2457 2289 2458 2459 2460 2290 // renew bad list 2291 @l = init_bads( MaxDeg ); 2292 @l = zAddBad( @l, @z, MaxDeg, 0 ); 2461 2293 2462 2463 2294 } else 2295 { 2464 2296 2465 2466 2467 @l = zAddBad( @l, @q, MaxDeg, 0 ); /// add all possible 'leadexps' !2468 2297 if ( size(@q) > 0 ) 2298 { 2299 @l = zAddBad( @l, @q, MaxDeg, 0 ); // add all possible 'leadexps' ! 2300 }; 2469 2301 2470 2471 2472 2473 2474 2475 2476 @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); /// kill from @tt all entries with leadexp in @l[@d]2477 2478 2479 2480 2481 2482 2483 2484 2485 }; 2486 2487 / //////////////////////////////////////////////////////////////////////////////2302 }; 2303 2304 PBW[index(k)] = PBW_next( PBW, k, SYS ); 2305 2306 PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k1)] , @l ); 2307 2308 @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); // kill from @tt all entries with leadexp in @l[@d] 2309 2310 k ++; 2311 }; 2312 2313 QNF_done(); 2314 my_var_done(); 2315 2316 return (makeIdeal(@z)); 2317 }; 2318 2319 /******************************************************/ 2488 2320 static proc center_vectorspace( int MaxDeg ) 2489 "2321 " 2490 2322 pure calculation of center as a finitely dimmensional Vector Space (deg <= MaxDeg ) 2491 2323 " 2492 2324 { 2493 2494 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 2505 2506 }; 2507 2508 / //////////////////////////////////////////////////////////////////////////////2325 my_var_init(); 2326 2327 int N = nvars( basering ); 2328 list P = PBW_base( MaxDeg, my_var(1) ); 2329 2330 for( int v = 2; v <= N; v++ ) 2331 { 2332 P = calc_k_base( P, my_var(v) ); 2333 }; 2334 2335 my_var_done(); 2336 2337 return( calc_k_base ( P ) ); 2338 }; 2339 2340 /******************************************************/ 2509 2341 /* static */ proc center_min_vectorspace( int MaxDeg ) 2510 "2342 " 2511 2343 computes the 'minimal' set of central elements (of deg <= MaxDeg) 2512 2344 by reducing the set of it's generators as vector space … … 2517 2349 { 2518 2350 2519 2520 2521 2522 2523 2524 }; 2525 2526 2527 / //////////////////////////////////////////////////////////////////////////////2528 / //////////////////////////////////////////////////////////////////////////////2529 // /centralizer's computations2530 2531 / //////////////////////////////////////////////////////////////////////////////2351 QNF_init(); 2352 ideal res = makeIdeal( zReduce( center_vectorspace( MaxDeg))); 2353 QNF_done(); 2354 2355 return( res ); 2356 }; 2357 2358 2359 /******************************************************/ 2360 /******************************************************/ 2361 // centralizer's computations 2362 2363 /******************************************************/ 2532 2364 static proc centralizer_vectorspace( poly p, int MaxDeg ) 2533 2365 { 2534 2535 }; 2536 2537 / //////////////////////////////////////////////////////////////////////////////2366 return( calc_k_base( PBW_base( MaxDeg, p ))); 2367 }; 2368 2369 /******************************************************/ 2538 2370 /* static */ proc centralizer_min_vectorspace( poly p, int MaxDeg ) 2539 2371 { 2540 2372 2541 2542 2543 2544 2545 2546 }; 2547 2548 / //////////////////////////////////////////////////////////////////////////////2373 QNF_init(); 2374 ideal res = makeIdeal( zReduce( centralizer_vectorspace( p, MaxDeg))); 2375 QNF_done(); 2376 2377 return( res ); 2378 }; 2379 2380 /******************************************************/ 2549 2381 /* static */ proc centralizer_min_iterative( poly p, int MaxDeg, list # ) 2550 "2382 " 2551 2383 computes the 'minimal' set of elements (of deg <= MaxDeg) generating centralizer of p in iterative way 2552 2384 Note: based on calc_k_base … … 2556 2388 " 2557 2389 { 2558 QNF_init(); 2559 2560 int n = myInt(#); 2561 int m = (MaxDeg < 0); 2562 2563 int MinDeg = 6; /// starting guess for MaxDeg 2564 int Delta = 4; /// increment of MaxDeg 2565 2566 if( m ) 2567 { 2568 // minimal guess 2569 MaxDeg = MinDeg; 2570 }; 2571 2572 list @q; 2573 2574 list PBW = list(); 2575 list PBW_PLAIN = list(); 2576 2577 PBW[ index(0) ] = PBW_init(); 2578 2579 list SYS = PBW_sys( my_vars(), p ); 2580 2581 list @z = list (); // result list 2582 list @l = init_bads( MaxDeg ); // verbotten loeadexps... 2583 2584 @q = PBW[ index(0) ]; 2585 2586 for (int k = 1; k <= ( MaxDeg+1 ); k ++ ) 2587 { 2588 @q = zRefine( calc_k_base(@q), 1 ); 2589 2590 if ( size(@q) > 0 ) 2591 { 2592 @z = @z + @q; /// computed desired elements 2593 2594 if( (n > 0) and (size(@z) > n) ) 2595 { 2596 break; // found needed elements 2597 }; 2598 }; 2599 2600 if( k == ( MaxDeg+1 ) ) 2601 { 2602 if( (n == 0) or ( !m ) ) 2603 { 2604 break; /// that's all 2605 }; 2606 2607 MaxDeg = MaxDeg + Delta; 2608 2609 // renew bad list 2610 @l = init_bads( MaxDeg ); 2611 @l = zAddBad( @l, @z, MaxDeg, 0 ); 2612 2613 } else 2614 { 2390 QNF_init(); 2391 2392 int n = myInt(#); 2393 int m = (MaxDeg < 0); 2394 2395 int MinDeg = 6; // starting guess for MaxDeg 2396 int Delta = 4; // increment of MaxDeg 2397 2398 if( m ) 2399 { 2400 // minimal guess 2401 MaxDeg = MinDeg; 2402 }; 2403 2404 list @q; 2405 2406 list PBW = list(); 2407 list PBW_PLAIN = list(); 2408 2409 PBW[ index(0) ] = PBW_init(); 2410 2411 list SYS = PBW_sys( my_vars(), p ); 2412 2413 list @z = list (); // result list 2414 list @l = init_bads( MaxDeg ); // verbotten loeadexps... 2415 2416 @q = PBW[ index(0) ]; 2417 2418 for (int k = 1; k <= ( MaxDeg+1 ); k ++ ) 2419 { 2420 @q = zRefine( calc_k_base(@q), 1 ); 2421 2422 if ( size(@q) > 0 ) 2423 { 2424 @z = @z + @q; // computed desired elements 2425 2426 if( (n > 0) and (size(@z) > n) ) 2427 { 2428 break; // found needed elements 2429 }; 2430 }; 2431 2432 if( k == ( MaxDeg+1 ) ) 2433 { 2434 if( (n == 0) or ( !m ) ) 2435 { 2436 break; // that's all 2437 }; 2438 2439 MaxDeg = MaxDeg + Delta; 2440 2441 // renew bad list 2442 @l = init_bads( MaxDeg ); 2443 @l = zAddBad( @l, @z, MaxDeg, 0 ); 2444 2445 } else 2446 { 2447 2448 if ( size(@q) > 0 ) 2449 { 2450 @l = zAddBad( @l, @q, MaxDeg, 0 ); // add all possible 'leadexps' ! 2451 }; 2452 2453 }; 2615 2454 2616 if ( size(@q) > 0 ) 2617 { 2618 @l = zAddBad( @l, @q, MaxDeg, 0 ); /// add all possible 'leadexps' ! 2619 }; 2620 2621 }; 2622 2623 PBW[index(k)] = PBW_next( PBW, k, SYS ); 2624 PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k1)] , @l ); 2625 @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); // kill from @tt all entries with leadexp in @l[@d] 2626 2627 }; 2628 2629 QNF_done(); 2630 return (makeIdeal(@z)); 2631 }; 2632 2633 2634 2635 /////////////////////////////////////////////////////////////////////////////// 2636 /////////////////////////////////////////////////////////////////////////////// 2455 PBW[index(k)] = PBW_next( PBW, k, SYS ); 2456 PBW_PLAIN = PBW_PLAIN + zCancel( PBW[index(k1)] , @l ); 2457 @q = PBW_PLAIN + zCancel( PBW[index(k)] , @l ); // kill from @tt all entries with leadexp in @l[@d] 2458 2459 }; 2460 2461 QNF_done(); 2462 return (makeIdeal(@z)); 2463 }; 2464 2465 2466 2467 /******************************************************/ 2468 /******************************************************/ 2637 2469 // main aliases 2638 2470 2639 /****************************************************** ************************/2471 /******************************************************/ 2640 2472 proc inCenter( def a ) 2641 "USAGE: inCenter(a); a poly/list/ideal2473 "USAGE: inCenter(a); a poly/list/ideal 2642 2474 RETURN: integer (1 if a in center, 0 otherwise) 2643 2475 EXAMPLE: example inCenter; shows examples" 2644 2476 { 2645 2646 2647 2648 2649 2650 2651 2652 2653 2654 2655 2656 2657 2658 2659 2660 2661 2662 2663 2477 QNF_init(); 2478 def res; 2479 2480 if ( typeof(a) == "poly" ) 2481 { 2482 res = inCenter_poly(a); 2483 } else 2484 { 2485 if ( (typeof(a)=="list") or (typeof(a)=="ideal") ) 2486 { 2487 res = inCenter_list(a); 2488 } else 2489 { 2490 res = a; 2491 }; 2492 }; 2493 2494 QNF_done(); 2495 return (res); 2664 2496 } 2665 2497 example … … 2685 2517 /******************************************************************************/ 2686 2518 proc inCentralizer( def a, poly f ) 2687 "USAGE: inCentralizer(a, f); a poly/list/ideal, f poly2519 "USAGE: inCentralizer(a, f); a poly/list/ideal, f poly 2688 2520 RETURN: integer (1 if a in centralizer(f), 0 otherwise) 2689 2521 EXAMPLE: example inCentralizer; shows examples" 2690 2522 { 2691 2692 2693 2694 2695 2696 2697 2698 2699 2700 2701 2702 2703 2704 2705 2706 2707 2708 2709 2523 QNF_init(); 2524 def res; 2525 2526 if ( typeof(a) == "poly" ) 2527 { 2528 res = inCentralizer_poly(a, f); 2529 } else 2530 { 2531 if ( (typeof(a)=="list") or (typeof(a)=="ideal") ) 2532 { 2533 res = inCentralizer_list(a, f); 2534 } else 2535 { 2536 res = a; 2537 }; 2538 }; 2539 2540 QNF_done(); 2541 return (res); 2710 2542 } 2711 2543 example … … 2719 2551 ncalgebra(1,D); // this is U(sl_2) 2720 2552 poly f = x^2; 2721 poly a = z; // /lies in center2553 poly a = z; // lies in center 2722 2554 poly b = y^2; 2723 2555 inCentralizer(a, f); … … 2731 2563 2732 2564 2733 / //////////////////////////////////////////////////////////////////////////////2565 /******************************************************/ 2734 2566 proc center(int MaxDeg, list # ) 2735 "USAGE: center(MaxDeg[, N]); int MaxDeg, int N2567 "USAGE: center(MaxDeg[, N]); int MaxDeg, int N 2736 2568 RETURN: ideal generated by found elements 2737 2569 NOTE: computes the 'minimal' set of central elements. … … 2744 2576 { 2745 2577 2746 if ( myInt(#) > 0 ) /// given a number of central elements to compute2747 2748 2749 2750 2751 if( MaxDeg >= 0 ) /// standard computation  the fastest one (often)2752 2753 //return ( center_min_iterative( MaxDeg, # ) );2754 2755 2756 2757 2758 2578 if ( myInt(#) > 0 ) // given a number of central elements to compute 2579 { 2580 return ( center_min_iterative( MaxDeg, # ) ); 2581 }; 2582 2583 if( MaxDeg >= 0 ) // standard computation  the fastest one (often) 2584 { 2585 // return ( center_min_iterative( MaxDeg, # ) ); 2586 return ( center_min_vectorspace ( MaxDeg ) ); 2587 }; 2588 2589 print( "Error: wrong arguments." ); 2590 return(); 2759 2591 2760 2592 } 2761 2593 example 2762 2594 { "EXAMPLE:"; echo = 2; 2763 2764 2765 2766 2767 2768 2769 ideal Z = center(2); /// find all central elements of degree <= 22770 2771 2772 ideal ZZ = center(1, 1 ); /// find the first non trivial central element2773 ZZ;2774 2775 }; 2776 2777 2778 / //////////////////////////////////////////////////////////////////////////////2595 ring a=0,(x,y,z),dp; 2596 matrix D[3][3]=0; 2597 D[1,2]=z; 2598 D[1,3]=2*x; 2599 D[2,3]=2*y; 2600 ncalgebra(1,D); // this is U(sl_2) 2601 ideal Z = center(2); // find all central elements of degree <= 2 2602 Z; 2603 inCenter(Z); 2604 ideal ZZ = center(1, 1 ); // find the first non trivial central element 2605 ZZ; ""; 2606 inCenter(ZZ); 2607 }; 2608 2609 2610 /******************************************************/ 2779 2611 proc centralizer( poly p, int MaxDeg, list # ) 2780 "USAGE: centralizer(F, MaxDeg[, N]); poly F, int MaxDeg, int N2612 "USAGE: centralizer(F, MaxDeg[, N]); poly F, int MaxDeg, int N 2781 2613 RETURN: ideal generated by found elements 2782 2614 NOTE: computes the 'minimal' set of elements of centralizer(F). … … 2789 2621 { 2790 2622 2791 2792 2793 2794 2795 2796 2797 2798 //return( centralizer_min_iterative( p, MaxDeg, # ) );2799 2800 2801 2802 2803 2623 if( myInt(#) > 0 ) 2624 { 2625 return( centralizer_min_iterative( p, MaxDeg, # ) ); 2626 }; 2627 2628 if( MaxDeg >= 0 ) 2629 { 2630 // return( centralizer_min_iterative( p, MaxDeg, # ) ); 2631 return( centralizer_min_vectorspace( p, MaxDeg ) ); 2632 }; 2633 2634 print( "Error: wrong arguments." ); 2635 return(); 2804 2636 2805 2637 } 2806 2638 example 2807 2639 { "EXAMPLE:"; echo = 2; 2808 ring a=0,(x,y,z),dp; 2809 matrix D[3][3]=0; 2810 D[1,2]=z; 2811 D[1,3]=2*x; 2812 D[2,3]=2*y; 2813 ncalgebra(1,D); // this is U(sl_2) 2814 poly f = 4*x*y+z^22*z; /// central polynomial 2815 f; 2816 ideal c = centralizer(f, 2); /// find all elements of degree <= 2 which lies in centralizer of f 2817 c; 2818 inCentralizer(c, f); 2819 ideal cc = centralizer(f, 1, 2 ); /// find at least first two non trivial elements of the centralizer of f 2820 cc; 2821 inCentralizer(cc, f); 2822 poly g = z^22*z; /// any polynomial 2823 g; 2824 c = centralizer(g, 2); /// find all elements of degree <= 2 which lies in centralizer of g 2825 c; 2826 inCentralizer(c, g); 2827 cc = centralizer(g, 1, 2 ); /// find at least first two non trivial elements of the centralizer of g 2828 cc; 2829 inCentralizer(cc, g); 2830 }; 2831 2832 2833 /////////////////////////////////////////////////////////////////////////////// 2834 /////////////////////////////////////////////////////////////////////////////// 2835 /////////////////////////////////////////////////////////////////////////////// 2640 ring a=0,(x,y,z),dp; 2641 matrix D[3][3]=0; 2642 D[1,2]=z; 2643 D[1,3]=2*x; 2644 D[2,3]=2*y; 2645 ncalgebra(1,D); // this is U(sl_2) 2646 poly f = 4*x*y+z^22*z; // central polynomial 2647 f; 2648 ideal c = centralizer(f, 2); // find all elements of degree <= 2 which lies in centralizer of f 2649 c; 2650 inCentralizer(c, f); 2651 ideal cc = centralizer(f, 1, 2 ); // find at least first two non trivial elements of the centralizer of f 2652 cc; 2653 inCentralizer(cc, f); 2654 poly g = z^22*z; // any polynomial 2655 g; ""; 2656 c = centralizer(g, 2); // find all elements of degree <= 2 which lies in centralizer of g 2657 c; ""; 2658 inCentralizer(c, g); 2659 cc = centralizer(g, 1, 2 ); // find at least first two non trivial elements of the centralizer of g 2660 cc; ""; 2661 inCentralizer(cc, g); 2662 };
Note: See TracChangeset
for help on using the changeset viewer.