Changeset a2c0302 in git
- Timestamp:
- Apr 22, 1997, 5:40:32 PM (27 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 56f7c59395cec760c41b7311172e225a5cc099eb
- Parents:
- eef3aacee99ae0cf329b82837ba229a249b16451
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
factory/fac_univar.cc
reef3aa ra2c0302 1 1 // emacs edit mode for this file is -*- C++ -*- 2 // $Id: fac_univar.cc,v 1. 6 1997-04-08 10:33:19schmidt Exp $2 // $Id: fac_univar.cc,v 1.7 1997-04-22 15:40:32 schmidt Exp $ 3 3 4 4 /* 5 5 $Log: not supported by cvs2svn $ 6 Revision 1.6 1997/04/08 10:33:19 schmidt 7 #include <config.h> added 8 6 9 Revision 1.5 1997/03/27 09:54:41 schmidt 7 10 timing output changed to TIMING … … 30 33 31 34 */ 35 32 36 33 37 #include <config.h> … … 53 57 TIMING_DEFINE_PRINT(fac_combineFactors); 54 58 55 #define MAX_FP_FAC 3 59 60 const int max_fp_fac = 3; 56 61 57 62 static modpk theModulus; 58 63 64 // !!! this should be placed in cf_gcd.h 59 65 CanonicalForm 60 66 iextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b ); 61 67 68 62 69 #ifdef DEBUGOUTPUT 70 #define DEBOUTHPRINT(stream, msg, hg) \ 71 {stream << deb_level_msg << msg, stream.flush(); hprint( hg ); stream << endl;} 72 63 73 static void 64 74 hprint ( int * a ) … … 72 82 i++; 73 83 } 74 cerr << ")" << endl << endl; 75 } 76 #endif 84 cerr << ")"; 85 } 86 #else /* DEBUGOUTPUT */ 87 #define DEBOUTHPRINT(stream, msg, hg) 88 #endif /* DEBUGOUTPUT */ 77 89 78 90 static void … … 123 135 for ( i = F; i.hasItem(); i++ ) 124 136 if ( (k = i.getItem().factor().degree()) < n ) 125 if ( k != 0 ) 137 if ( k == -1 ) { 138 WARN( k == -1, "there occured an error. factory was not able to factorize\n" 139 "correctly mod p. Please send the example which caused\n" 140 "this error to the authors. Nonetheless we will go on with the\n" 141 "calculations hoping the result will be correct. Thank you." ); 142 } 143 else if ( k != 0 ) 126 144 a[k] = 1; 127 145 } … … 134 152 for ( i = 1; i < m; i++ ) 135 153 if ( (k = F[i].degree()) < n ) 136 if ( k != 0 ) 154 if ( k == -1 ) { 155 WARN( k == -1, "there occured an error. factory was not able to factorize\n" 156 "correctly mod p. Please send the example which caused\n" 157 "this error to the authors. Nonetheless we will go on with the\n" 158 "calculations hoping the result will be correct. Thank you." ); 159 } 160 else if ( k != 0 ) 137 161 a[k] = 1; 138 162 } … … 167 191 CanonicalForm sum = 0; 168 192 for ( i = f; i.hasTerms(); i++ ) sum += i.coeff() * i.coeff(); 169 DEBOUTLN( cerr, "sum = " ,sum );193 DEBOUTLN( cerr, "sum = " << sum ); 170 194 return sqrt( cf2double( sum ) ); 171 195 } … … 174 198 kBound ( const CanonicalForm & f, int p ) 175 199 { 176 DEBOUTLN( cerr, "lc(f) = " ,lc(f) );177 return (int)((f.degree()*log(2)+log( fabs(cf2double(lc(f))) )+log( dnorm( f ) )) / log( p ) + 0.5) + 1;200 DEBOUTLN( cerr, "lc(f) = " << lc(f) ); 201 return (int)((f.degree()*log(2)+log( fabs(cf2double(lc(f))) )+log( dnorm( f ) )) / log( (double)p ) + 0.5) + 1; 178 202 } 179 203 … … 190 214 return false; 191 215 else if ( degree( f ) + degree( theFactors[i] ) == d ) { 192 DEBOUTLN( cerr, "ldfr (f) = " ,f );193 DEBOUTLN( cerr, "ldfr (g) = " ,theFactors[i] );216 DEBOUTLN( cerr, "ldfr (f) = " << f ); 217 DEBOUTLN( cerr, "ldfr (g) = " << theFactors[i] ); 194 218 CanonicalForm g = pp( pk( recip_lf * f * theFactors[i] ) ); 195 DEBOUTLN( cerr, "ldfr (pk(f*g)) = " ,g );219 DEBOUTLN( cerr, "ldfr (pk(f*g)) = " << g ); 196 220 CanonicalForm gg, hh; 197 DEBOUTLN( cerr, "F = " ,F );198 DEBOUTLN( cerr, "g = " ,g );221 DEBOUTLN( cerr, "F = " << F ); 222 DEBOUTLN( cerr, "g = " << g ); 199 223 if ( divremt( F, g, gg, hh ) && hh.isZero() ) { 200 224 ZF.append( CFFactor( g, exp ) ); … … 219 243 } 220 244 245 221 246 static int 222 247 choosePrimes ( int * p, const CanonicalForm & f ) … … 227 252 int prime; 228 253 229 while ( ptr < maxp && i < MAX_FP_FAC) {254 while ( ptr < maxp && i < max_fp_fac ) { 230 255 prime = cf_getPrime( ptr ); 231 256 if ( mod( lc( f ), prime ) != 0 ) { … … 239 264 ptr++; 240 265 } 241 return ( i == MAX_FP_FAC);266 return ( i == max_fp_fac ); 242 267 } 243 268 … … 250 275 int i, j, save; 251 276 int p = pk.getp(), k = pk.getk(); 252 int no_iter = (int)(log( k)/log(2)+2);277 int no_iter = (int)(log( (double)k )/log(2)+2); 253 278 int * kvals = new int[no_iter]; 254 279 255 DEBOUTSL( cerr ); 256 DEBOUT( cerr, "quadratic lift called with p = " << p << " and k = " << k ); 257 DEBOUTENDL( cerr ); 280 DEBOUTLN( cerr, "quadratic lift called with p = " << p << " and k = " << k ); 258 281 for ( j = 0, i = k; i > 1; i = ( i+1 ) / 2, j++ ) kvals[j] = i; 259 282 kvals[j] = 1; … … 282 305 j--; 283 306 setCharacteristic( p, kvals[j+1] ); 284 DEBOUTSL( cerr ); 285 DEBOUT( cerr, "lifting from p^" << kvals[j+1] << " to p^" << kvals[j] ); 286 DEBOUTENDL( cerr ); 307 DEBOUTLN( cerr, "lifting from p^" << kvals[j+1] << " to p^" << kvals[j] ); 287 308 c = mapinto( c ); 288 DEBOUTLN( cerr, " !!! g = " ,mapinto( g ) );309 DEBOUTLN( cerr, " !!! g = " << mapinto( g ) ); 289 310 g1 = mapinto( lf ) / mapinto( lc( g ) ) * mapinto( g ); 290 311 h1 = mapinto( lf ) / mapinto( lc( h ) ) * mapinto( h ); 291 312 // (void)extgcd( g1, h1, a, b ); 292 // DEBOUTLN( cerr, " a = " ,aa );293 // DEBOUTLN( cerr, " b = " ,bb );313 // DEBOUTLN( cerr, " a = " << aa ); 314 // DEBOUTLN( cerr, " b = " << bb ); 294 315 a = mapinto( a ); b = mapinto( b ); 295 316 a += ( ( 1 - a * g1 ) * a ) % h1; 296 317 b += ( ( 1 - b * h1 ) * b ) % g1; 297 DEBOUTLN( cerr, " a = " ,a );298 DEBOUTLN( cerr, " b = " ,b );318 DEBOUTLN( cerr, " a = " << a ); 319 DEBOUTLN( cerr, " b = " << b ); 299 320 divrem( a * c, h1, q, r ); 300 321 tmp = b * c + q * g1; … … 308 329 modulus = power( CanonicalForm(p), kvals[j] ); 309 330 if ( mod( f - g * h, modulus ) != 0 ) 310 DEBOUTLN( cerr, "error at lift stage " ,i );331 DEBOUTLN( cerr, "error at lift stage " << i ); 311 332 i++; 312 333 } … … 375 396 } 376 397 398 377 399 CFFList 378 400 ZFactorizeUnivariate( const CanonicalForm& ff, bool issqrfree ) … … 383 405 int i, k, exp, n; 384 406 bool ok; 385 CFFList H, G, F[MAX_FP_FAC];407 CFFList H, F[max_fp_fac]; 386 408 CFFList ZF; 387 int * p = new int [ MAX_FP_FAC];409 int * p = new int [max_fp_fac]; 388 410 int * D = 0; 389 411 int * Dh = 0; 390 412 ListIterator<CFFactor> J, I; 413 414 DEBINCLEVEL( cerr, "ZFactorizeUnivariate" ); 391 415 On( SW_SYMMETRIC_FF ); 416 417 // get squarefree decomposition of f 392 418 if ( issqrfree ) 393 419 H.append( CFFactor( g, 1 ) ); 394 420 else 395 421 H = sqrFree( g ); 422 423 DEBOUTLN( cerr, "H = " << H ); 424 425 // cycle through squarefree factors of f 396 426 for ( J = H; J.hasItem(); ++J ) { 397 427 f = J.getItem().factor(); 398 428 if ( f.inCoeffDomain() ) continue; 399 429 n = f.degree() / 2 + 1; 400 if ( D != 0 ) { 401 delete [] D; 402 delete [] Dh; 403 } 430 delete [] D; 431 delete [] Dh; 404 432 D = new int [n]; D[0] = n; 405 433 Dh = new int [n]; Dh[0] = n; 406 434 exp = J.getItem().exp(); 435 436 // choose primes to factor f 407 437 TIMING_START(fac_choosePrimes); 408 438 ok = choosePrimes( p, f ); 409 439 TIMING_END_AND_PRINT(fac_choosePrimes, "time to choose the primes: "); 410 440 if ( ! ok ) { 411 DEBOUTLN( cerr, "error: no good prime found to factorize ", f ); 412 ASSERT( 0, "error: no good prime found" ); 441 DEBOUTLN( cerr, "warning: no good prime found to factorize " << f ); 442 WARN( ok, "there occured an error. We went out of primes p\n" 443 "to factorize mod p. Please send the example which caused\n" 444 "this error to the authors. Nonetheless we will go on with the\n" 445 "calculations hoping the result will be correct. Thank you."); 413 446 ZF.append( CFFactor( f, exp ) ); 414 } 415 else { 416 TIMING_START(fac_facModPrimes); 417 for ( i = 0; i < MAX_FP_FAC; i++ ) { 418 setCharacteristic( p[i] ); 419 fp = mapinto( f ); 420 F[i] = FpFactorizeUnivariateCZ( fp, true ); 447 continue; 448 } 449 450 // factorize f modulo certain primes 451 TIMING_START(fac_facModPrimes); 452 for ( i = 0; i < max_fp_fac; i++ ) { 453 setCharacteristic( p[i] ); 454 fp = mapinto( f ); 455 F[i] = FpFactorizeUnivariateCZ( fp, true ); 421 456 // if ( p[i] < 23 && fp.degree() < 10 ) 422 457 // F[i] = FpFactorizeUnivariateB( fp, true ); 423 458 // else 424 459 // F[i] = FpFactorizeUnivariateCZ( fp, true ); 425 DEBOUTSL( cerr ); 426 DEBOUT( cerr, "F[i] = " << F[i] << ", p = " << p[i] ); 427 DEBOUTENDL( cerr ); 460 DEBOUTLN( cerr, "F[i] = " << F[i] << ", p = " << p[i] ); 461 } 462 TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: "); 463 setCharacteristic( 0 ); 464 465 // do some strange things with the D's 466 initHG( D, F[0] ); 467 hgroup( D ); 468 DEBOUTHPRINT( cerr, "D = ", D ); 469 for ( i = 1; i < max_fp_fac; i++ ) { 470 initHG( Dh, F[i] ); 471 hgroup( Dh ); 472 DEBOUTHPRINT( cerr, "Dh = ", Dh ); 473 hintersect( D, Dh ); 474 DEBOUTHPRINT( cerr, "D = ", D ); 475 } 476 477 // look which p gives the shortest factorization of f mod p 478 // j: index of that p in p[] 479 int min, j; 480 min = F[0].length(), j = 0; 481 for ( i = 1; i < max_fp_fac; i++ ) { 482 if ( min >= F[i].length() ) { 483 j = i; min = F[i].length(); 428 484 } 429 TIMING_END_AND_PRINT(fac_facModPrimes, "time to factorize mod primes: "); 430 setCharacteristic( 0 ); 431 #ifdef DEBUGOUTPUT 432 DEBOUTLN( cerr, "D = ", ' ' ); 433 hprint( D ); 434 #endif 435 initHG( D, F[0] ); 436 hgroup( D ); 437 #ifdef DEBUGOUTPUT 438 DEBOUTLN( cerr, "D = ", ' ' ); 439 hprint( D ); 440 #endif 441 for ( i = 1; i < MAX_FP_FAC; i++ ) { 442 initHG( Dh, F[i] ); 443 hgroup( Dh ); 444 #ifdef DEBUGOUTPUT 445 DEBOUTLN( cerr, "Dh = ", ' ' ); 446 hprint( Dh ); 447 #endif 448 hintersect( D, Dh ); 449 #ifdef DEBUGOUTPUT 450 DEBOUTLN( cerr, "D = ", ' ' ); 451 hprint( D ); 452 #endif 453 454 } 455 int min, j; 456 min = F[0].length(), j = 0; 457 for ( i = 1; i < MAX_FP_FAC; i++ ) { 458 if ( min >= F[i].length() ) { 459 j = i; min = F[i].length(); 460 } 461 } 462 k = kBound( f, p[j] ); 463 CFArray theFactors( F[j].length() ); 485 } 486 k = kBound( f, p[j] ); 487 CFArray theFactors( F[j].length() ); 464 488 // pk = power( CanonicalForm( p[j] ), k ); 465 489 // pkhalf = pk / 2; 466 467 DEBOUTLN( cerr, "coeff bound = ",pk.getpk() );468 469 470 471 472 473 474 475 DEBOUTLN( cerr, "factor to lift = ",I.getItem().factor() );476 477 478 479 480 481 482 DEBOUTLN( cerr, "dummy1 = ",dummy1 );483 DEBOUTLN( cerr, "dummy2 = ",dummy2 );484 485 486 487 488 489 490 DEBOUTLN( cerr, "F[j] = ",F[j] );491 492 493 DEBOUTLN( cerr, "i = ",i );494 DEBOUTLN( cerr, "dummy1 = ",dummy1 );495 490 modpk pk( p[j], k ); 491 DEBOUTLN( cerr, "coeff bound = " << pk.getpk() ); 492 theModulus = pk; 493 setCharacteristic( p[j] ); 494 fp = mapinto( f ); 495 F[j].sort( cmpFactor ); 496 I = F[j]; i = 0; 497 TIMING_START(fac_liftFactors); 498 while ( I.hasItem() ) { 499 DEBOUTLN( cerr, "factor to lift = " << I.getItem().factor() ); 500 if ( isOn( SW_FAC_QUADRATICLIFT ) ) 501 ok = UnivariateQuadraticLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 ); 502 else 503 ok = UnivariateLinearLift( f, I.getItem().factor(), fp / I.getItem().factor(), pk, lc( f ), dummy1, dummy2 ); 504 if ( ok ) { 505 // should be done in a more efficient way 506 DEBOUTLN( cerr, "dummy1 = " << dummy1 ); 507 DEBOUTLN( cerr, "dummy2 = " << dummy2 ); 508 f = dummy2; 509 fp /= I.getItem().factor(); 510 ZF.append( CFFactor( dummy1, exp ) ); 511 I.remove( 0 ); 512 I = F[j]; 513 i = 0; 514 DEBOUTLN( cerr, "F[j] = " << F[j] ); 515 } 516 else { 517 DEBOUTLN( cerr, "i = " << i ); 518 DEBOUTLN( cerr, "dummy1 = " << dummy1 ); 519 setCharacteristic( 0 ); 496 520 // theFactors[i] = pk( dummy1 * pk.inverse( lc( dummy1 ) ) ); 497 theFactors[i] = pk( dummy1 ); 498 setCharacteristic( p[j] ); 499 i++; 500 I++; 501 } 521 theFactors[i] = pk( dummy1 ); 522 setCharacteristic( p[j] ); 523 i++; 524 I++; 502 525 } 503 TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: "); 504 DEBOUTLN( cerr, "ZF = ", ZF ); 505 initHG( Dh, theFactors ); 506 hgroup( Dh ); 507 #ifdef DEBUGOUTPUT 508 DEBOUTLN( cerr, "Dh = ", ' ' ); 509 hprint( Dh ); 510 #endif 511 hintersect( D, Dh ); 512 setCharacteristic( 0 ); 513 for ( int l = i; l < F[j].length(); l++ ) 514 theFactors[l] = 1; 515 DEBOUTLN( cerr, "theFactors = ", theFactors ); 516 DEBOUTLN( cerr, "f = ", f ); 517 DEBOUTSL( cerr ); 518 DEBOUT( cerr, "p = " << pk.getp() << ", k = " << pk.getk() ); 519 DEBOUTENDL( cerr ); 520 #ifdef DEBUGOUTPUT 521 DEBOUTLN( cerr, "D = ", ' ' ); 522 hprint( D ); 523 #endif 524 lf = lc( f ); 525 (void)iextgcd( pk.getpk(), lf, dummy1, recip_lf ); 526 DEBOUTLN( cerr, "recip_lf = ", recip_lf ); 527 TIMING_START(fac_combineFactors); 528 for ( i = 1; i < D[0]; i++ ) 529 if ( D[i] != 0 ) 530 while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) ); 531 if ( degree( f ) > 0 ) 532 ZF.append( CFFactor( f, exp ) ); 533 TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: "); 534 } 535 } 526 } 527 TIMING_END_AND_PRINT(fac_liftFactors, "time to lift the factors: "); 528 DEBOUTLN( cerr, "ZF = " << ZF ); 529 initHG( Dh, theFactors ); 530 hgroup( Dh ); 531 DEBOUTHPRINT( cerr, "Dh = ", Dh ); 532 hintersect( D, Dh ); 533 setCharacteristic( 0 ); 534 for ( int l = i; l < F[j].length(); l++ ) 535 theFactors[l] = 1; 536 DEBOUTLN( cerr, "theFactors = " << theFactors ); 537 DEBOUTLN( cerr, "f = " << f ); 538 DEBOUTLN( cerr, "p = " << pk.getp() << ", k = " << pk.getk() ); 539 DEBOUTHPRINT( cerr, "D = ", D ); 540 lf = lc( f ); 541 (void)iextgcd( pk.getpk(), lf, dummy1, recip_lf ); 542 DEBOUTLN( cerr, "recip_lf = " << recip_lf ); 543 TIMING_START(fac_combineFactors); 544 for ( i = 1; i < D[0]; i++ ) 545 if ( D[i] != 0 ) 546 while ( liftDegreeFactRec( theFactors, f, recip_lf, lf, pk, 0, i, ZF, exp ) ); 547 if ( degree( f ) > 0 ) 548 ZF.append( CFFactor( f, exp ) ); 549 TIMING_END_AND_PRINT(fac_combineFactors, "time to combine the factors: "); 550 } 551 552 // brush up our result 536 553 if ( ZF.getFirst().factor().inCoeffDomain() ) 537 554 ZF.removeFirst(); … … 540 557 else 541 558 ZF.insert( CFFactor( cont, 1 ) ); 542 if ( D != 0 ) { 543 delete [] D; 544 delete [] Dh; 545 } 559 delete [] D; 560 delete [] Dh; 546 561 if ( ! symmsave ) 547 562 Off( SW_SYMMETRIC_FF ); 563 564 DEBDECLEVEL( cerr, "ZFactorizeUnivariate" ); 548 565 return ZF; 549 566 }
Note: See TracChangeset
for help on using the changeset viewer.