Changeset dd2a9b3 in git
 Timestamp:
 May 13, 2011, 9:44:24 AM (12 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0604212ebb110535022efecad887940825b97c3f')
 Children:
 a2dd9b278a186a85e5955f8e7adad56483fddffd
 Parents:
 4775bdf78134ad29d5eb334ec54f2e1fa7f460ff
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

factory/facHensel.cc
r4775bd rdd2a9b3 204 204 } 205 205 206 zz_pX kronSubFp (const CanonicalForm& A, int d) 207 { 208 int degAy= degree (A); 209 zz_pX result; 210 result.rep.SetLength (d*(degAy + 1)); 211 212 zz_p *resultp; 213 resultp= result.rep.elts(); 214 zz_pX buf; 215 zz_p *bufp; 216 217 for (CFIterator i= A; i.hasTerms(); i++) 218 { 219 if (i.coeff().inCoeffDomain()) 220 buf= convertFacCF2NTLzzpX (i.coeff()); 221 else 222 buf= convertFacCF2NTLzzpX (i.coeff()); 223 224 int k= i.exp()*d; 225 bufp= buf.rep.elts(); 226 int bufRepLength= (int) buf.rep.length(); 227 for (int j= 0; j < bufRepLength; j++) 228 resultp [j + k]= bufp [j]; 229 } 230 result.normalize(); 231 232 return result; 233 } 234 235 zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha) 236 { 237 int degAy= degree (A); 238 zz_pEX result; 239 result.rep.SetLength (d*(degAy + 1)); 240 241 Variable v; 242 zz_pE *resultp; 243 resultp= result.rep.elts(); 244 zz_pEX buf1; 245 zz_pE *buf1p; 246 zz_pX buf2; 247 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha)); 248 249 for (CFIterator i= A; i.hasTerms(); i++) 250 { 251 if (i.coeff().inCoeffDomain()) 252 { 253 buf2= convertFacCF2NTLzzpX (i.coeff()); 254 buf1= to_zz_pEX (to_zz_pE (buf2)); 255 } 256 else 257 buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo); 258 259 int k= i.exp()*d; 260 buf1p= buf1.rep.elts(); 261 int buf1RepLength= (int) buf1.rep.length(); 262 for (int j= 0; j < buf1RepLength; j++) 263 resultp [j + k]= buf1p [j]; 264 } 265 result.normalize(); 266 267 return result; 268 } 269 270 void 271 kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d, 272 const Variable& alpha) 273 { 274 int degAy= degree (A); 275 subA1.rep.SetLength ((long) d*(degAy + 1)); 276 subA2.rep.SetLength ((long) d*(degAy + 1)); 277 278 Variable v; 279 zz_pE *subA1p; 280 zz_pE *subA2p; 281 subA1p= subA1.rep.elts(); 282 subA2p= subA2.rep.elts(); 283 zz_pEX buf; 284 zz_pE *bufp; 285 zz_pX buf2; 286 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha)); 287 288 for (CFIterator i= A; i.hasTerms(); i++) 289 { 290 if (i.coeff().inCoeffDomain()) 291 { 292 buf2= convertFacCF2NTLzzpX (i.coeff()); 293 buf= to_zz_pEX (to_zz_pE (buf2)); 294 } 295 else 296 buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo); 297 298 int k= i.exp()*d; 299 int kk= (degAy  i.exp())*d; 300 bufp= buf.rep.elts(); 301 int bufRepLength= (int) buf.rep.length(); 302 for (int j= 0; j < bufRepLength; j++) 303 { 304 subA1p [j + k]= bufp [j]; 305 subA2p [j + kk]= bufp [j]; 306 } 307 } 308 subA1.normalize(); 309 subA2.normalize(); 310 } 311 312 void 313 kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d) 314 { 315 int degAy= degree (A); 316 subA1.rep.SetLength ((long) d*(degAy + 1)); 317 subA2.rep.SetLength ((long) d*(degAy + 1)); 318 319 Variable v; 320 zz_p *subA1p; 321 zz_p *subA2p; 322 subA1p= subA1.rep.elts(); 323 subA2p= subA2.rep.elts(); 324 zz_pX buf; 325 zz_p *bufp; 326 327 for (CFIterator i= A; i.hasTerms(); i++) 328 { 329 buf= convertFacCF2NTLzzpX (i.coeff()); 330 331 int k= i.exp()*d; 332 int kk= (degAy  i.exp())*d; 333 bufp= buf.rep.elts(); 334 int bufRepLength= (int) buf.rep.length(); 335 for (int j= 0; j < bufRepLength; j++) 336 { 337 subA1p [j + k]= bufp [j]; 338 subA2p [j + kk]= bufp [j]; 339 } 340 } 341 subA1.normalize(); 342 subA2.normalize(); 343 } 344 345 CanonicalForm 346 reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k, 347 const Variable& alpha) 348 { 349 Variable y= Variable (2); 350 Variable x= Variable (1); 351 352 zz_pEX f= F; 353 zz_pEX g= G; 354 int degf= deg(f); 355 int degg= deg(g); 356 357 zz_pEX buf1; 358 zz_pEX buf2; 359 zz_pEX buf3; 360 361 zz_pE *buf1p; 362 zz_pE *buf2p; 363 zz_pE *buf3p; 364 if (f.rep.length() < (long) d*(k+1)) //zero padding 365 f.rep.SetLength ((long)d*(k+1)); 366 367 zz_pE *gp= g.rep.elts(); 368 zz_pE *fp= f.rep.elts(); 369 CanonicalForm result= 0; 370 int i= 0; 371 int lf= 0; 372 int lg= d*k; 373 int degfSubLf= degf; 374 int deggSubLg= degglg; 375 int repLengthBuf2; 376 int repLengthBuf1; 377 int ii; 378 zz_pE zzpEZero= zz_pE(); 379 380 while (degf >= lf  lg >= 0) 381 { 382 if (degfSubLf >= d) 383 repLengthBuf1= d; 384 else if (degfSubLf < 0) 385 repLengthBuf1= 0; 386 else 387 repLengthBuf1= degfSubLf + 1; 388 buf1.rep.SetLength((long) repLengthBuf1); 389 390 buf1p= buf1.rep.elts(); 391 for (int ind= 0; ind < repLengthBuf1; ind++) 392 buf1p [ind]= fp [ind + lf]; 393 buf1.normalize(); 394 395 repLengthBuf1= buf1.rep.length(); 396 397 if (deggSubLg >= d  1) 398 repLengthBuf2= d  1; 399 else if (deggSubLg < 0) 400 repLengthBuf2= 0; 401 else 402 repLengthBuf2= deggSubLg + 1; 403 404 buf2.rep.SetLength ((long) repLengthBuf2); 405 buf2p= buf2.rep.elts(); 406 for (int ind= 0; ind < repLengthBuf2; ind++) 407 { 408 buf2p [ind]= gp [ind + lg]; 409 } 410 buf2.normalize(); 411 412 repLengthBuf2= buf2.rep.length(); 413 414 buf3.rep.SetLength((long) repLengthBuf2 + d); 415 buf3p= buf3.rep.elts(); 416 buf2p= buf2.rep.elts(); 417 buf1p= buf1.rep.elts(); 418 for (int ind= 0; ind < repLengthBuf1; ind++) 419 buf3p [ind]= buf1p [ind]; 420 for (int ind= repLengthBuf1; ind < d; ind++) 421 buf3p [ind]= zzpEZero; 422 for (int ind= 0; ind < repLengthBuf2; ind++) 423 buf3p [ind + d]= buf2p [ind]; 424 buf3.normalize(); 425 426 result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i); 427 i++; 428 429 430 lf= i*d; 431 degfSubLf= degf  lf; 432 433 lg= d*(ki); 434 deggSubLg= degg  lg; 435 436 buf1p= buf1.rep.elts(); 437 438 if (lg >= 0 && deggSubLg > 0) 439 { 440 if (repLengthBuf2 > degfSubLf + 1) 441 degfSubLf= repLengthBuf2  1; 442 int tmp= tmin (repLengthBuf1, deggSubLg); 443 for (int ind= 0; ind < tmp; ind++) 444 gp [ind + lg] = buf1p [ind]; 445 } 446 447 if (lg < 0) 448 break; 449 450 buf2p= buf2.rep.elts(); 451 if (degfSubLf >= 0) 452 { 453 for (int ind= 0; ind < repLengthBuf2; ind++) 454 fp [ind + lf] = buf2p [ind]; 455 } 456 } 457 458 return result; 459 } 460 461 CanonicalForm 462 reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k) 463 { 464 Variable y= Variable (2); 465 Variable x= Variable (1); 466 467 zz_pX f= F; 468 zz_pX g= G; 469 int degf= deg(f); 470 int degg= deg(g); 471 472 zz_pX buf1; 473 zz_pX buf2; 474 zz_pX buf3; 475 476 zz_p *buf1p; 477 zz_p *buf2p; 478 zz_p *buf3p; 479 480 if (f.rep.length() < (long) d*(k+1)) //zero padding 481 f.rep.SetLength ((long)d*(k+1)); 482 483 zz_p *gp= g.rep.elts(); 484 zz_p *fp= f.rep.elts(); 485 CanonicalForm result= 0; 486 int i= 0; 487 int lf= 0; 488 int lg= d*k; 489 int degfSubLf= degf; 490 int deggSubLg= degglg; 491 int repLengthBuf2; 492 int repLengthBuf1; 493 int ii; 494 zz_p zzpZero= zz_p(); 495 while (degf >= lf  lg >= 0) 496 { 497 if (degfSubLf >= d) 498 repLengthBuf1= d; 499 else if (degfSubLf < 0) 500 repLengthBuf1= 0; 501 else 502 repLengthBuf1= degfSubLf + 1; 503 buf1.rep.SetLength((long) repLengthBuf1); 504 505 buf1p= buf1.rep.elts(); 506 for (int ind= 0; ind < repLengthBuf1; ind++) 507 buf1p [ind]= fp [ind + lf]; 508 buf1.normalize(); 509 510 repLengthBuf1= buf1.rep.length(); 511 512 if (deggSubLg >= d  1) 513 repLengthBuf2= d  1; 514 else if (deggSubLg < 0) 515 repLengthBuf2= 0; 516 else 517 repLengthBuf2= deggSubLg + 1; 518 519 buf2.rep.SetLength ((long) repLengthBuf2); 520 buf2p= buf2.rep.elts(); 521 for (int ind= 0; ind < repLengthBuf2; ind++) 522 buf2p [ind]= gp [ind + lg]; 523 524 buf2.normalize(); 525 526 repLengthBuf2= buf2.rep.length(); 527 528 529 buf3.rep.SetLength((long) repLengthBuf2 + d); 530 buf3p= buf3.rep.elts(); 531 buf2p= buf2.rep.elts(); 532 buf1p= buf1.rep.elts(); 533 for (int ind= 0; ind < repLengthBuf1; ind++) 534 buf3p [ind]= buf1p [ind]; 535 for (int ind= repLengthBuf1; ind < d; ind++) 536 buf3p [ind]= zzpZero; 537 for (int ind= 0; ind < repLengthBuf2; ind++) 538 buf3p [ind + d]= buf2p [ind]; 539 buf3.normalize(); 540 541 result += convertNTLzzpX2CF (buf3, x)*power (y, i); 542 i++; 543 544 545 lf= i*d; 546 degfSubLf= degf  lf; 547 548 lg= d*(ki); 549 deggSubLg= degg  lg; 550 551 buf1p= buf1.rep.elts(); 552 553 if (lg >= 0 && deggSubLg > 0) 554 { 555 if (repLengthBuf2 > degfSubLf + 1) 556 degfSubLf= repLengthBuf2  1; 557 int tmp= tmin (repLengthBuf1, deggSubLg); 558 for (int ind= 0; ind < tmp; ind++) 559 gp [ind + lg] = buf1p [ind]; 560 } 561 if (lg < 0) 562 break; 563 564 buf2p= buf2.rep.elts(); 565 if (degfSubLf >= 0) 566 { 567 for (int ind= 0; ind < repLengthBuf2; ind++) 568 fp [ind + lf] = buf2p [ind]; 569 } 570 } 571 572 return result; 573 } 574 575 CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha) 576 { 577 Variable y= Variable (2); 578 Variable x= Variable (1); 579 580 zz_pEX f= F; 581 zz_pE *fp= f.rep.elts(); 582 583 zz_pEX buf; 584 zz_pE *bufp; 585 CanonicalForm result= 0; 586 int i= 0; 587 int degf= deg(f); 588 int k= 0; 589 int degfSubK; 590 int repLength; 591 while (degf >= k) 592 { 593 degfSubK= degf  k; 594 if (degfSubK >= d) 595 repLength= d; 596 else 597 repLength= degfSubK + 1; 598 599 buf.rep.SetLength ((long) repLength); 600 bufp= buf.rep.elts(); 601 for (int j= 0; j < repLength; j++) 602 bufp [j]= fp [j + k]; 603 buf.normalize(); 604 605 result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i); 606 i++; 607 k= d*i; 608 } 609 610 return result; 611 } 612 613 CanonicalForm reverseSubstFp (const zz_pX& F, int d) 614 { 615 Variable y= Variable (2); 616 Variable x= Variable (1); 617 618 zz_pX f= F; 619 zz_p *fp= f.rep.elts(); 620 621 zz_pX buf; 622 zz_p *bufp; 623 CanonicalForm result= 0; 624 int i= 0; 625 int degf= deg(f); 626 int k= 0; 627 int degfSubK; 628 int repLength; 629 while (degf >= k) 630 { 631 degfSubK= degf  k; 632 if (degfSubK >= d) 633 repLength= d; 634 else 635 repLength= degfSubK + 1; 636 637 buf.rep.SetLength ((long) repLength); 638 bufp= buf.rep.elts(); 639 for (int j= 0; j < repLength; j++) 640 bufp [j]= fp [j + k]; 641 buf.normalize(); 642 643 result += convertNTLzzpX2CF (buf, x)*power (y, i); 644 i++; 645 k= d*i; 646 } 647 648 return result; 649 } 650 651 // assumes input to be reduced mod M and to be an element of Fq not Fp 652 CanonicalForm 653 mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const 654 CanonicalForm& M) 655 { 656 int d1= tmax (degree (F, 1), degree (G, 1)) + 1; 657 int d2= tmax (degree (F, 2), degree (G, 2)); 658 659 zz_pX F1, F2; 660 kronSubRecipro (F1, F2, F, d1); 661 zz_pX G1, G2; 662 kronSubRecipro (G1, G2, G, d1); 663 664 int k= d1*degree (M); 665 MulTrunc (F1, F1, G1, (long) k); 666 667 mul (F2, F2, G2); 668 F2 >>= k; 669 670 return reverseSubst (F1, F2, d1, d2); 671 } 672 673 //Kronecker substitution 674 CanonicalForm 675 mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const 676 CanonicalForm& M) 677 { 678 CanonicalForm A= F; 679 CanonicalForm B= G; 680 681 int p= getCharacteristic(); 682 683 int degAx= degree (A, 1); 684 int degAy= degree (A, 2); 685 int degBx= degree (B, 1); 686 int degBy= degree (B, 2); 687 int d1= degAx + 1 + degBx; 688 int d2= tmax (degree (A, 2), degree (B, 2)); 689 690 if (d1 > 128 && d2 > 160 && (degAy == degBy)) 691 return mulMod2NTLFpReci (A, B, M); 692 693 zz_pX NTLA= kronSubFp (A, d1); 694 zz_pX NTLB= kronSubFp (B, d1); 695 696 int k= d1*degree (M); 697 MulTrunc (NTLA, NTLA, NTLB, (long) k); 698 699 A= reverseSubstFp (NTLA, d1); 700 701 return A; 702 } 703 704 // assumes input to be reduced mod M and to be an element of Fq not Fp 705 CanonicalForm 706 mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const 707 CanonicalForm& M, const Variable& alpha) 708 { 709 int d1= tmax (degree (F, 1), degree (G, 1)) + 1; 710 int d2= tmax (degree (F, 2), degree (G, 2)); 711 712 zz_pEX F1, F2; 713 kronSubRecipro (F1, F2, F, d1, alpha); 714 zz_pEX G1, G2; 715 kronSubRecipro (G1, G2, G, d1, alpha); 716 717 int k1= d1*degree (M); 718 MulTrunc (F1, F1, G1, (long) k1); 719 720 mul (F2, F2, G2); 721 722 F2 >>= k1; 723 724 CanonicalForm result= reverseSubst (F1, F2, d1, d2, alpha); 725 726 return result; 727 } 728 729 CanonicalForm 730 mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const 731 CanonicalForm& M) 732 { 733 Variable alpha; 734 CanonicalForm A= F; 735 CanonicalForm B= G; 736 737 if (hasFirstAlgVar (A, alpha)  hasFirstAlgVar (B, alpha)) 738 { 739 int degAx= degree (A, 1); 740 int degAy= degree (A, 2); 741 int degBx= degree (B, 1); 742 int degBy= degree (B, 2); 743 int d1= degAx + degBx + 1; 744 int d2= tmax (degree (A, 2), degree (B, 2)); 745 zz_p::init (getCharacteristic()); 746 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha)); 747 zz_pE::init (NTLMipo); 748 749 int degMipo= degree (getMipo (alpha)); 750 if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy)) 751 return mulMod2NTLFqReci (A, B, M, alpha); 752 753 zz_pEX NTLA= kronSub (A, d1, alpha); 754 zz_pEX NTLB= kronSub (B, d1, alpha); 755 756 int k= d1*degree (M); 757 758 MulTrunc (NTLA, NTLA, NTLB, (long) k); 759 760 A= reverseSubst (NTLA, d1, alpha); 761 762 return A; 763 } 764 else 765 return mulMod2NTLFp (A, B, M); 766 } 767 206 768 CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B, 207 769 const CanonicalForm& M) … … 220 782 int degG= degree (G, y); 221 783 222 if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) && 784 if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) && 223 785 (F.level() == G.level())) 224 { 786 { 225 787 CanonicalForm result= mulNTL (F, G); 226 788 return mod (result, M); 227 789 } 228 790 else if (degF <= 1 && degG <= 1) 229 { 791 { 230 792 CanonicalForm result= F*G; 231 793 return mod (result, M); 232 794 } 795 796 int sizeF= size (F); 797 int sizeG= size (G); 798 799 int fallBackToNaive= 50; 800 if (sizeF < fallBackToNaive  sizeG < fallBackToNaive) 801 return mod (F*G, M); 802 803 if (CFFactory::gettype() != GaloisFieldDomain && 804 (((degFdegG) < 50 && degF > degG)  ((degGdegF) < 50 && degF <= degG))) 805 return mulMod2NTLFq (F, G, M); 233 806 234 807 int m= (int) ceil (degree (M)/2.0); … … 404 977 return F; 405 978 CanonicalForm A= F; 406 Variable y= F.mvar();979 Variable y= Variable (2); 407 980 Variable x= Variable (1); 408 A= swapvar (A, x, y); 409 CanonicalForm result= 0; 410 CFIterator i= A; 411 while (d  i.exp() < 0) 412 i++; 413 414 for (; i.hasTerms() && (d  i.exp() >= 0); i++) 415 result += swapvar (i.coeff(),x,y)*power (x, d  i.exp()); 416 return result; 981 if (degree (A, x) > 0) 982 { 983 A= swapvar (A, x, y); 984 CanonicalForm result= 0; 985 CFIterator i= A; 986 while (d  i.exp() < 0) 987 i++; 988 989 for (; i.hasTerms() && (d  i.exp() >= 0); i++) 990 result += swapvar (i.coeff(),x,y)*power (x, d  i.exp()); 991 return result; 992 } 993 else 994 return A*power (x, d); 417 995 } 418 996 … … 490 1068 491 1069 CanonicalForm Q; 492 if (degB <= 1 )1070 if (degB <= 1  CFFactory::gettype() == GaloisFieldDomain) 493 1071 { 494 1072 CanonicalForm R; … … 526 1104 } 527 1105 528 if (degB <= 1 )1106 if (degB <= 1  CFFactory::gettype() == GaloisFieldDomain) 529 1107 { 530 1108 divrem2 (A, B, Q, R, M);
Note: See TracChangeset
for help on using the changeset viewer.