Changeset c244d04 in git for factory/fac_ezgcd.cc
 Timestamp:
 Apr 24, 2012, 1:01:15 PM (11 years ago)
 Branches:
 (u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')
 Children:
 9f84ad2f3376850c4d78697edb203d32c07d3094
 Parents:
 f7a4e9271997208f4ee40a14fd1e1b5c6c0c1ffc
 gitauthor:
 Martin Lee <martinlee84@web.de>20120424 13:01:15+02:00
 gitcommitter:
 Martin Lee <martinlee84@web.de>20120507 14:15:53+02:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/fac_ezgcd.cc
rf7a4e9 rc244d04 1 /* emacs edit mode for this file is * C++ * */ 2 /* $Id$ */ 1 /*****************************************************************************\ 2 * Computer Algebra System SINGULAR 3 \*****************************************************************************/ 4 /** @file fac_ezgcd.cc 5 * 6 * This file implements the GCD of two multivariate polynomials over Q using 7 * EZGCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor, 8 * Labahnn 9 * 10 * @author Martin Lee 11 * 12 **/ 13 /*****************************************************************************/ 3 14 4 15 #include "config.h" … … 9 20 #include "cf_defs.h" 10 21 #include "canonicalform.h" 11 #include "fac_util.h" 22 #include "fac_util.h" // HEADER for this file 12 23 #include "cf_algorithm.h" 13 24 #include "cf_reval.h" … … 331 342 } 332 343 333 static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG ); 334 335 static CanonicalForm ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal ); 336 337 static CanonicalForm ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound ); 338 339 static modpk findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG ); 340 341 static modpk enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk ); 344 static CanonicalForm 345 ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, 346 bool internal ) 347 { 348 bool isRat= isOn (SW_RATIONAL); 349 if (!isRat) 350 On (SW_RATIONAL); 351 CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, 352 lcD, cand, result; 353 CFArray DD( 1, 2 ), lcDD( 1, 2 ); 354 int degF, degG, delta, t; 355 REvaluation bt; 356 bool gcdfound = false; 357 Variable x = Variable(1); 358 modpk bound; 359 360 F= FF; 361 G= GG; 362 363 CFMap M,N; 364 int smallestDegLev; 365 int best_level= compress4EZGCD (F, G, M, N, smallestDegLev); 366 367 if (best_level == 0) 368 { 369 if (!isRat) 370 Off (SW_RATIONAL); 371 return G.genOne(); 372 } 373 374 F= M (F); 375 G= M (G); 376 377 378 DEBINCLEVEL( cerr, "ezgcd" ); 379 DEBOUTLN( cerr, "FF = " << FF ); 380 DEBOUTLN( cerr, "GG = " << GG ); 381 f = content( F, x ); g = content( G, x ); d = gcd( f, g ); 382 DEBOUTLN( cerr, "f = " << f ); 383 DEBOUTLN( cerr, "g = " << g ); 384 F /= f; G /= g; 385 if ( F.isUnivariate() && G.isUnivariate() ) 386 { 387 DEBDECLEVEL( cerr, "ezgcd" ); 388 if (!isRat) 389 Off (SW_RATIONAL); 390 if(F.mvar()==G.mvar()) 391 d*=gcd(F,G); 392 return N (d); 393 } 394 else if ( gcd_test_one( F, G, false ) ) 395 { 396 DEBDECLEVEL( cerr, "ezgcd" ); 397 if (!isRat) 398 Off (SW_RATIONAL); 399 return N (d); 400 } 401 lcF = LC( F, x ); lcG = LC( G, x ); 402 lcD = gcd( lcF, lcG ); 403 delta = 0; 404 degF = degree( F, x ); degG = degree( G, x ); 405 t = tmax( F.level(), G.level() ); 406 //bound = findBound( F, G, lcF, lcG, degF, degG ); 407 if ( ! internal ) 408 b = REvaluation( 2, t, IntRandom( 25 ) ); 409 while ( ! gcdfound ) 410 { 411 /// > A2 412 DEBOUTLN( cerr, "search for evaluation, delta = " << delta ); 413 DEBOUTLN( cerr, "F = " << F ); 414 DEBOUTLN( cerr, "G = " << G ); 415 findeval( F, G, Fb, Gb, Db, b, delta, degF, degG ); 416 DEBOUTLN( cerr, "found evaluation b = " << b ); 417 DEBOUTLN( cerr, "F(b) = " << Fb ); 418 DEBOUTLN( cerr, "G(b) = " << Gb ); 419 DEBOUTLN( cerr, "D(b) = " << Db ); 420 delta = degree( Db ); 421 /// > A3 422 if ( delta == 0 ) 423 { 424 DEBDECLEVEL( cerr, "ezgcd" ); 425 if (!isRat) 426 Off (SW_RATIONAL); 427 return N (d); 428 } 429 /// > A4 430 //deltaold = delta; 431 while ( 1 ) 432 { 433 bt = b; 434 findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG ); 435 int dd=degree( Dbt ); 436 if ( dd /*degree( Dbt )*/ == 0 ) 437 { 438 DEBDECLEVEL( cerr, "ezgcd" ); 439 if (!isRat) 440 Off (SW_RATIONAL); 441 return N (d); 442 } 443 if ( dd /*degree( Dbt )*/ == delta ) 444 break; 445 else if ( dd /*degree( Dbt )*/ < delta ) 446 { 447 delta = dd /*degree( Dbt )*/; 448 b = bt; 449 Db = Dbt; Fb = Fbt; Gb = Gbt; 450 } 451 } 452 DEBOUTLN( cerr, "now after A4, delta = " << delta ); 453 /// > A5 454 if (delta == degF) 455 { 456 if (degF <= degG && fdivides (F, G)) 457 { 458 DEBDECLEVEL( cerr, "ezgcd" ); 459 if (!isRat) 460 Off (SW_RATIONAL); 461 return N (d*F); 462 } 463 else 464 delta; 465 } 466 else if (delta == degG) 467 { 468 if (degG <= degF && fdivides( G, F )) 469 { 470 DEBDECLEVEL( cerr, "ezgcd" ); 471 if (!isRat) 472 Off (SW_RATIONAL); 473 return N (d*G); 474 } 475 else 476 delta; 477 } 478 if ( delta != degF && delta != degG ) 479 { 480 /// > A6 481 bool B_is_F; 482 CanonicalForm xxx1, xxx2; 483 CanonicalForm buf; 484 DD[1] = Fb / Db; 485 buf= Gb/Db; 486 xxx1 = gcd( DD[1], Db ); 487 xxx2 = gcd( buf, Db ); 488 if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) && 489 (size (F) <= size (G))) 490  (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain())) 491 { 492 B = F; 493 DD[2] = Db; 494 lcDD[1] = lcF; 495 lcDD[2] = lcD; 496 B_is_F = true; 497 } 498 else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) && 499 (size (G) < size (F))) 500  (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain())) 501 { 502 DD[1] = buf; 503 B = G; 504 DD[2] = Db; 505 lcDD[1] = lcG; 506 lcDD[2] = lcD; 507 B_is_F = false; 508 } 509 else 510 { 511 //special case 512 Off(SW_USE_EZGCD); 513 result=gcd( F, G ); 514 On(SW_USE_EZGCD); 515 if (!isRat) 516 Off (SW_RATIONAL); 517 return N (d*result); 518 } 519 /// > A7 520 DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) ); 521 DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) ); 522 DEBOUTLN( cerr, "(hensel) B = " << B ); 523 DEBOUTLN( cerr, "(hensel) lcB = " << LC( B, Variable(1) ) ); 524 DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) ); 525 DEBOUTLN( cerr, "(hensel) DD = " << DD ); 526 DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD ); 527 gcdfound= Hensel (B*lcD, DD, b, lcDD); 528 DEBOUTLN( cerr, "(hensel finished) DD = " << DD ); 529 530 if (gcdfound == 1) 531 { 532 Off (SW_USE_EZGCD); 533 result= gcd (F,G); 534 On (SW_USE_EZGCD); 535 if (!isRat) 536 Off (SW_RATIONAL); 537 return N (d*result); 538 } 539 540 if (gcdfound) 541 { 542 cand = DD[2] / content( DD[2], Variable(1) ); 543 gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); 544 } 545 /// > A8 (gcdfound) 546 } 547 delta; 548 } 549 /// > A9 550 DEBDECLEVEL( cerr, "ezgcd" ); 551 if (!isRat) 552 Off (SW_RATIONAL); 553 return N (d*cand); 554 } 342 555 343 556 CanonicalForm 344 557 ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG ) 345 558 { 346 REvaluation b; 347 return ezgcd( FF, GG, b, false ); 348 } 349 350 static CanonicalForm 351 ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal ) 352 { 353 CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD, cand, result; 354 CFArray DD( 1, 2 ), lcDD( 1, 2 ); 355 int degF, degG, delta, t; 356 REvaluation bt; 357 bool gcdfound = false; 358 Variable x = Variable(1); 359 modpk bound; 360 361 F= FF; 362 G= GG; 363 364 CFMap M,N; 365 int smallestDegLev; 366 int best_level= compress4EZGCD (F, G, M, N, smallestDegLev); 367 368 if (best_level == 0) return G.genOne(); 369 370 F= M (F); 371 G= M (G); 372 373 374 DEBINCLEVEL( cerr, "ezgcd" ); 375 DEBOUTLN( cerr, "FF = " << FF ); 376 DEBOUTLN( cerr, "GG = " << GG ); 377 f = content( F, x ); g = content( G, x ); d = gcd( f, g ); 378 DEBOUTLN( cerr, "f = " << f ); 379 DEBOUTLN( cerr, "g = " << g ); 380 F /= f; G /= g; 381 if ( F.isUnivariate() && G.isUnivariate() ) 382 { 383 DEBDECLEVEL( cerr, "ezgcd" ); 384 if(F.mvar()==G.mvar()) 385 d*=gcd(F,G); 386 return N (d); 387 } 388 else if ( gcd_test_one( F, G, false ) ) 389 { 390 DEBDECLEVEL( cerr, "ezgcd" ); 391 return N (d); 392 } 393 lcF = LC( F, x ); lcG = LC( G, x ); 394 lcD = gcd( lcF, lcG ); 395 delta = 0; 396 degF = degree( F, x ); degG = degree( G, x ); 397 t = tmax( F.level(), G.level() ); 398 //bound = findBound( F, G, lcF, lcG, degF, degG ); 399 if ( ! internal ) 400 b = REvaluation( 2, t, IntRandom( 25 ) ); 401 while ( ! gcdfound ) 402 { 403 /// > A2 404 DEBOUTLN( cerr, "search for evaluation, delta = " << delta ); 405 DEBOUTLN( cerr, "F = " << F ); 406 DEBOUTLN( cerr, "G = " << G ); 407 findeval( F, G, Fb, Gb, Db, b, delta, degF, degG ); 408 DEBOUTLN( cerr, "found evaluation b = " << b ); 409 DEBOUTLN( cerr, "F(b) = " << Fb ); 410 DEBOUTLN( cerr, "G(b) = " << Gb ); 411 DEBOUTLN( cerr, "D(b) = " << Db ); 412 delta = degree( Db ); 413 /// > A3 414 if ( delta == 0 ) 415 { 416 DEBDECLEVEL( cerr, "ezgcd" ); 417 return N (d); 418 } 419 /// > A4 420 //deltaold = delta; 421 while ( 1 ) { 422 bt = b; 423 findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG ); 424 int dd=degree( Dbt ); 425 if ( dd /*degree( Dbt )*/ == 0 ) 426 { 427 DEBDECLEVEL( cerr, "ezgcd" ); 428 return N (d); 429 } 430 if ( dd /*degree( Dbt )*/ == delta ) 431 break; 432 else if ( dd /*degree( Dbt )*/ < delta ) { 433 delta = dd /*degree( Dbt )*/; 434 b = bt; 435 Db = Dbt; Fb = Fbt; Gb = Gbt; 436 } 437 } 438 DEBOUTLN( cerr, "now after A4, delta = " << delta ); 439 /// > A5 440 if (delta == degF) 441 { 442 if (degF <= degG && fdivides (F, G)) 443 { 444 DEBDECLEVEL( cerr, "ezgcd" ); 445 return N (d*F); 446 } 447 else 448 delta; 449 } 450 else if (delta == degG) 451 { 452 if (degG <= degF && fdivides( G, F )) 453 { 454 DEBDECLEVEL( cerr, "ezgcd" ); 455 return N (d*G); 456 } 457 else 458 delta; 459 } 460 if ( delta != degF && delta != degG ) 461 { 462 /// > A6 463 bool B_is_F; 464 CanonicalForm xxx1, xxx2; 465 CanonicalForm buf; 466 DD[1] = Fb / Db; 467 buf= Gb/Db; 468 xxx1 = gcd( DD[1], Db ); 469 xxx2 = gcd( buf, Db ); 470 if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) && 471 (size (F) <= size (G))) 472  (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain())) 473 { 474 B = F; 475 DD[2] = Db; 476 lcDD[1] = lcF; 477 lcDD[2] = lcD; 478 B_is_F = true; 479 } 480 else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) && 481 (size (G) < size (F))) 482  (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain())) 483 { 484 DD[1] = buf; 485 B = G; 486 DD[2] = Db; 487 lcDD[1] = lcG; 488 lcDD[2] = lcD; 489 B_is_F = false; 490 } 491 else 492 { 493 Off(SW_USE_EZGCD); 494 result=gcd( F, G ); 495 On(SW_USE_EZGCD); 496 return N (d*result); 497 } 498 /*CanonicalForm xxx; 499 //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) { 500 DD[1] = Fb / Db; 501 xxx= gcd( DD[1], Db ); 502 DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx ); 503 DEBOUTLN( cerr, "Fb/Db = " << DD[1] ); 504 DEBOUTLN( cerr, "Db = " << Db ); 505 if (xxx.inCoeffDomain()) { 506 B = F; 507 DD[2] = Db; 508 lcDD[1] = lcF; 509 lcDD[2] = lcF; 510 B *= lcF; 511 B_is_F=true; 512 } 513 //else if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) { 514 else 515 { 516 DD[1] = Gb / Db; 517 xxx=gcd( DD[1], Db ); 518 DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx ); 519 DEBOUTLN( cerr, "Gb/Db = " << DD[1] ); 520 DEBOUTLN( cerr, "Db = " << Db ); 521 if (xxx.inCoeffDomain()) 522 { 523 B = G; 524 DD[2] = Db; 525 lcDD[1] = lcG; 526 lcDD[2] = lcG; 527 B *= lcG; 528 B_is_F=false; 529 } 530 else 531 { 532 #ifdef DEBUGOUTPUT 533 CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound ); 534 DEBDECLEVEL( cerr, "ezgcd" ); 535 return dummyres; 536 #else 537 return d * ezgcd_specialcase( F, G, b, bound ); 538 #endif 539 } 540 }*/ 541 /// > A7 542 DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) ); 543 DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) ); 544 //bound = enlargeBound( B, DD[2], DD[1], bound ); 545 DEBOUTLN( cerr, "(hensel) B = " << B ); 546 DEBOUTLN( cerr, "(hensel) lcB = " << LC( B, Variable(1) ) ); 547 DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) ); 548 DEBOUTLN( cerr, "(hensel) DD = " << DD ); 549 DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD ); 550 gcdfound= Hensel (B*lcD, DD, b, lcDD); 551 DEBOUTLN( cerr, "(hensel finished) DD = " << DD ); 552 553 if (gcdfound == 1) 554 { 555 Off (SW_USE_EZGCD); 556 result= gcd (F,G); 557 On (SW_USE_EZGCD_P); 558 return N (d*result); 559 } 560 561 if (gcdfound) 562 { 563 cand = DD[2] / content( DD[2], Variable(1) ); 564 gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); 565 } 566 /// > A8 (gcdfound) 567 } 568 delta; 569 } 570 /// > A9 571 DEBDECLEVEL( cerr, "ezgcd" ); 572 return N (d*cand); 573 } 574 575 static CanonicalForm 576 ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & /*b*/, const modpk & /*bound*/ ) 577 { 578 CanonicalForm d; 579 #if 1 580 Off(SW_USE_EZGCD); 581 //bool ntl0=isOn(SW_USE_NTL_GCD_0); 582 //Off(SW_USE_NTL_GCD_0); 583 //bool ntlp=isOn(SW_USE_NTL_GCD_P); 584 //Off(SW_USE_NTL_GCD_P); 585 d=gcd( F, G ); 586 //if (ntl0) On(SW_USE_NTL_GCD_0); 587 //if (ntlp) On(SW_USE_NTL_GCD_P); 588 On(SW_USE_EZGCD); 589 return d; 590 #else 591 /* 592 * I am not sure, if these lines are much of use. 593 * Remove? 594 */ 595 DEBOUTLN( cerr, "ezgcd: special case" ); 596 CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt; 597 CFArray DD( 1, 2 ), lcDD( 1, 2 ); 598 Variable x = Variable( 1 ); 599 bool gcdfound; 600 601 Dt = 1; 602 /// > S1 603 DEBOUTLN( cerr, "ezgcdspec: (S1)" ); 604 Ft = ezgcd( F, F.deriv( x ) ); 605 if ( Ft.isOne() ) 606 { 607 // In this case F is squarefree and we came here by bad chance 608 // (means: bad evaluation point). Try again with another 609 // evaluation point. Bug fix (?) by JS. The bad example was 610 // gcd.debug ocr /+USE_EZGCD/@12/CB 611 // '(16*B^8208*B^6*C+927*B^4*C^21512*B^2*C^3+432*C^4)' 612 // '(4*B^7*C^250*B^5*C^3+208*B^3*C^4288*B*C^5)' 613 b.nextpoint(); 614 return ezgcd( F, G, b, true ); 615 } 616 617 DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft ); 618 L = F / Ft; 619 /// > S2 620 DEBOUTLN( cerr, "ezgcdspec: (S2)" ); 621 622 L = ezgcd( L, G, b, true ); 623 DEBOUTLN( cerr, "ezgcdspec: (S2) done, Ds = " << Ds ); 624 Gt = G / L; 625 Ds = L; Dt = L; 626 Lb = b( L ); Gb = b( Gt ); Fb = b( Ft ); 627 d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) ); 628 while ( true ) { 629 /// > S3 630 DEBOUTLN( cerr, "ezgcdspec: (S3)" ); 631 DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb ); 632 DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb ); 633 DD[1] = gcd( Lb, gcd( Fb, Gb ) ); 634 /// > S4 635 DEBOUTLN( cerr, "ezgcdspec: (S4)" ); 636 if ( degree( DD[1] ) == 0 ) 637 return Ds; 638 /// > S5 639 DEBOUTLN( cerr, "ezgcdspec: (S5)" ); 640 DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb ); 641 DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] ); 642 Db = DD[1]; 643 if ( ! (DD[2] = Lb/DD[1]).isOne() ) { 644 DEBOUTLN( cerr, "ezgcdspec: (S55)" ); 645 lcDD[2] = LC( L, x ); 646 lcDD[1] = d; 647 DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] ); 648 DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] ); 649 LL = L * d; 650 modpk newbound = enlargeBound( LL, DD[2], DD[1], bound ); 651 DEBOUTLN( cerr, "ezgcdspec: begin with lift." ); 652 DEBOUTLN( cerr, "ezgcdspec: B = " << LL ); 653 DEBOUTLN( cerr, "ezgcdspec: DD = " << DD ); 654 DEBOUTLN( cerr, "ezgcdspec: lcDD = " << lcDD ); 655 DEBOUTLN( cerr, "ezgcdspec: b = " << b ); 656 DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() ); 657 DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) ); 658 DEBOUTLN( cerr, "ezgcdspec: test = " << b( LL )  DD[1] * DD[2] ); 659 gcdfound = Hensel( LL, DD, lcDD, b, newbound, x ); 660 ASSERT( gcdfound, "fatal error in algorithm" ); 661 Dt = pp( DD[1] ); 662 } 663 DEBOUTLN( cerr, "ezgcdspec: (S7)" ); 664 Ds = Ds * Dt; 665 Fb = Fb / Db; 666 Gb = Gb / Db; 667 } 668 // this point is never reached, but to avoid compiler warnings let's return a value 669 return 0; 670 #endif 671 } 672 673 static void 674 findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG ) 675 { 676 bool ok; 677 if ( delta != 0 ) 678 b.nextpoint(); 679 DEBOUTLN( cerr, "ezgcd: (findeval) F = " << F <<", G="<< G); 680 DEBOUTLN( cerr, "ezgcd: (findeval) degF = " << degF << ", degG="<<degG ); 681 do { 682 DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b ); 683 Fb = b( F ); 684 ok = degree( Fb ) == degF; 685 if ( ok ) { 686 Gb = b( G ); 687 ok = degree( Gb ) == degG; 688 } 689 690 if ( ok ) 691 { 692 Db = gcd( Fb, Gb ); 693 if ( delta > 0 ) 694 ok = degree( Db ) < delta; 695 } 696 if ( ! ok ) 697 { 698 b.nextpoint(); 699 } 700 } while ( ! ok ); 701 } 702 703 static modpk 704 enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk ) 705 { 706 DEBOUTLN( cerr, "ezgcd: (enlarge bound) F = " << F ); 707 DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb = " << Lb ); 708 DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db = " << Db ); 709 DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb % p = " << mod( Lb, pk.getp() ) ); 710 DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db % p = " << mod( Db, pk.getp() ) ); 711 712 CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) * 713 tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) ); 714 int p = pk.getp(); 715 int k = pk.getk(); 716 CanonicalForm bound = pk.getpk(); 717 while ( bound < limit ) { 718 k++; 719 bound *= p; 720 } 721 k *= 2; 722 DEBOUTLN( cerr, "ezgcd: (enlarge bound) newbound = " << power( CanonicalForm( p ), k ) ); 723 return modpk( p, k ); 724 } 725 726 static modpk 727 findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG ) 728 { 729 CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) * 730 gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) ); 731 int p, i = 0; 732 do { 733 p = cf_getBigPrime( i ); 734 i++; 735 } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() ); 736 CanonicalForm bound = p; 737 i = 1; 738 while ( bound < limit ) { 739 i++; 740 bound *= p; 741 } 742 return modpk( p, i ); 743 } 559 REvaluation b; 560 return ezgcd( FF, GG, b, false ); 561 } 562
Note: See TracChangeset
for help on using the changeset viewer.