Changeset 4e2cc1 in git
 Timestamp:
 Jul 11, 2014, 1:40:48 PM (9 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 03f640f7f62621edbfe49fd44dd3646e7f33da1c
 Parents:
 1a4c80b0c5ac3d9994f051b93db0759bf02ceb60
 gitauthor:
 Martin Lee <martinlee84@web.de>20140711 13:40:48+02:00
 gitcommitter:
 Martin Lee <martinlee84@web.de>20140717 12:07:11+02:00
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/cfModResultant.cc
r1a4c80b r4e2cc1 24 24 #include "cf_algorithm.h" 25 25 #include "cf_iter.h" 26 #include "cf_irred.h" 27 #include "cf_generator.h" 28 #include "cf_random.h" 29 #include "cf_map_ext.h" 30 #include "facFqBivarUtil.h" 26 31 27 32 #ifdef HAVE_NTL … … 253 258 if (F.isZero()  G.isZero()) 254 259 return 0; 260 Variable alpha; 255 261 256 262 #ifdef HAVE_FLINT 257 nmod_poly_t FLINTF, FLINTG; 258 convertFacCF2nmod_poly_t (FLINTF, F); 259 convertFacCF2nmod_poly_t (FLINTG, G); 260 mp_limb_t FLINTresult= nmod_poly_resultant (FLINTF, FLINTG); 261 nmod_poly_clear (FLINTF); 262 nmod_poly_clear (FLINTG); 263 return CanonicalForm ((long) FLINTresult); 263 if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha)) 264 { 265 nmod_poly_t FLINTF, FLINTG; 266 convertFacCF2nmod_poly_t (FLINTF, F); 267 convertFacCF2nmod_poly_t (FLINTG, G); 268 mp_limb_t FLINTresult= nmod_poly_resultant (FLINTF, FLINTG); 269 nmod_poly_clear (FLINTF); 270 nmod_poly_clear (FLINTG); 271 return CanonicalForm ((long) FLINTresult); 272 } 264 273 #else 274 if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha)) 275 { 276 if (fac_NTL_char != getCharacteristic()) 277 { 278 fac_NTL_char= getCharacteristic(); 279 zz_p::init (getCharacteristic()); 280 } 281 zz_pX NTLF= convertFacCF2NTLzzpX (F); 282 zz_pX NTLG= convertFacCF2NTLzzpX (G); 283 284 zz_p NTLResult= resultant (NTLF, NTLG); 285 286 return CanonicalForm (to_long (rep (NTLResult))); 287 } 288 #endif 289 //at this point F or G has an algebraic var. 265 290 if (fac_NTL_char != getCharacteristic()) 266 291 { … … 268 293 zz_p::init (getCharacteristic()); 269 294 } 270 zz_pX NTL F= convertFacCF2NTLzzpX (F);271 zz_p X NTLG= convertFacCF2NTLzzpX (G);272 273 zz_p NTLResult= resultant (NTLF, NTLG);274 275 return CanonicalForm (to_long (rep (NTLResult))); 276 #endif 295 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha)); 296 zz_pE::init (NTLMipo); 297 zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo); 298 zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo); 299 zz_pE NTLResult= resultant (NTLF, NTLG); 300 301 return convertNTLzzpE2CF (NTLResult, alpha); 277 302 #else 278 303 return resultant (F, G, F.mvar()); … … 282 307 static inline 283 308 void evalPoint (const CanonicalForm& F, const CanonicalForm& G, 284 CanonicalForm& FEval, CanonicalForm& GEval, int& evalPoint) 309 CanonicalForm& FEval, CanonicalForm& GEval, 310 CFGenerator& evalPoint) 285 311 { 286 312 int degF, degG; … … 290 316 do 291 317 { 292 evalPoint++; 293 if (evalPoint >= getCharacteristic()) 318 if (!evalPoint.hasItems()) 294 319 break; 295 FEval= F (evalPoint , 2);296 GEval= G (evalPoint , 2);320 FEval= F (evalPoint.item(), 2); 321 GEval= G (evalPoint.item(), 2); 297 322 if (degree (FEval, 1) < degF  degree (GEval, 1) < degG) 323 { 324 evalPoint.next(); 298 325 continue; 326 } 299 327 else 300 328 return; 301 329 } 302 while (evalPoint < getCharacteristic());330 while (evalPoint.hasItems()); 303 331 } 304 332 … … 350 378 Variable y= Variable (2); 351 379 352 int i= 1;353 380 CanonicalForm GEval, FEval, recResult, H; 354 381 CanonicalForm newtonPoly= 1; … … 358 385 int bound= degAx*degree (G, 2) + degree (F, 2)*degBx; 359 386 387 int p= getCharacteristic(); 388 CanonicalForm minpoly; 389 Variable alpha= Variable (tmax (F.level(), G.level()) + 1); 390 bool algExt= hasFirstAlgVar (F, alpha)  hasFirstAlgVar (G, alpha); 391 CFGenerator * gen; 392 bool extOfExt= false; 393 Variable v= alpha; 394 CanonicalForm primElemAlpha, imPrimElemAlpha; 395 CFList source,dest; 396 if (!algExt && (p < (1 << 28))) 397 { 398 // pass to an extension of size at least 2^29 399 // for very very large input that is maybe too small though 400 int deg= ceil (29.0*((double) log (2)/log (p)))+1; 401 minpoly= randomIrredpoly (deg, z); 402 alpha= rootOf (minpoly); 403 AlgExtGenerator AlgExtGen (alpha); 404 gen= AlgExtGen.clone(); 405 for (int i= 0; i < p; i++) // skip values from the prime field 406 (*gen).next(); 407 } 408 else if (!algExt) 409 { 410 FFGenerator FFGen; 411 gen= FFGen.clone(); 412 } 413 else 414 { 415 int deg= ceil (29.0*((double) log (2)/log (p))); 416 if (degree (getMipo (alpha)) < deg) 417 { 418 mpz_t field_size; 419 mpz_init (field_size); 420 mpz_ui_pow_ui (field_size, p, 421 deg + degree (getMipo (alpha))  deg%degree (getMipo (alpha))); 422 423 // field_size needs to fit in an int because of mapUp, mapDown, length of lists etc. 424 if (mpz_fits_sint_p (field_size)) 425 { 426 minpoly= randomIrredpoly (deg + degree (getMipo (alpha)) 427  deg%degree (getMipo (alpha)), z); 428 v= rootOf (minpoly); 429 Variable V_buf2; 430 bool primFail= false; 431 extOfExt= true; 432 primElemAlpha= primitiveElement (alpha, V_buf2, primFail); 433 ASSERT (!primFail, "failure in integer factorizer"); 434 if (primFail) 435 ; //ERROR 436 else 437 imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v); 438 F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest); 439 G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest); 440 } 441 else 442 { 443 deg= deg  deg % degree (getMipo (alpha)); 444 mpz_ui_pow_ui (field_size, p, deg); 445 while (deg / degree (getMipo (alpha)) >= 2 && !mpz_fits_sint_p (field_size)) 446 { 447 deg = degree (getMipo (alpha)); 448 mpz_ui_pow_ui (field_size, p, deg); 449 } 450 if (deg != degree (getMipo (alpha))) 451 { 452 minpoly= randomIrredpoly (deg, z); 453 v= rootOf (minpoly); 454 Variable V_buf2; 455 bool primFail= false; 456 extOfExt= true; 457 primElemAlpha= primitiveElement (alpha, V_buf2, primFail); 458 ASSERT (!primFail, "failure in integer factorizer"); 459 if (primFail) 460 ; //ERROR 461 else 462 imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v); 463 F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest); 464 G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest); 465 } 466 } 467 mpz_clear (field_size); 468 } 469 AlgExtGenerator AlgExtGen (v); 470 gen= AlgExtGen.clone(); 471 for (int i= 0; i < p; i++) 472 (*gen).next(); 473 } 360 474 int count= 0; 361 475 int equalCount= 0; 476 CanonicalForm point; 362 477 do 363 478 { 364 evalPoint (F, G, FEval, GEval, i); 365 366 ASSERT (i < getCharacteristic(), "ran out of points"); 367 368 recResult= resultantFp (FEval, GEval, z); 369 370 H= newtonInterp (i, recResult, newtonPoly, modResult, y); 479 evalPoint (F, G, FEval, GEval, *gen); 480 481 recResult= resultantFp (FEval, GEval, z, prob); 482 483 H= newtonInterp ((*gen).item(), recResult, newtonPoly, modResult, y); 371 484 372 485 if (H == modResult) … … 377 490 count++; 378 491 if (count > bound  (prob && equalCount == 2)) 379 break; 492 { 493 if (!algExt && degree (H, alpha) <= 0) 494 break; 495 else if (algExt) 496 { 497 if (extOfExt && !isInExtension (H, imPrimElemAlpha, 1, primElemAlpha, 498 dest, source)) 499 { 500 H= mapDown (H, primElemAlpha, imPrimElemAlpha, alpha, dest, source); 501 prune (v); 502 break; 503 } 504 else if (!extOfExt) 505 break; 506 } 507 } 380 508 381 509 modResult= H; 382 newtonPoly *= (y  i); 510 newtonPoly *= (y  (*gen).item()); 511 if ((*gen).hasItems()) 512 (*gen).next(); 513 else 514 STICKYASSERT (0, "out of evaluation points"); 383 515 } while (1); 516 517 delete gen; 384 518 385 519 return N (H);
Note: See TracChangeset
for help on using the changeset viewer.