1/************************************************************** 2 * 3 * factor.c 4 * 5 * General purpose factoring program 6 * 7 * Updates: 8 * 18 May 97 REC - invoked new, fast divide functions 9 * 26 Apr 97 RDW - fixed tabs and unix EOL 10 * 20 Apr 97 RDW - conversion to TC4.5 11 * 12 * c. 1997 Perfectly Scientific, Inc. 13 * All Rights Reserved. 14 * 15 * 16 *************************************************************/ 17 18/* include files */ 19 20#include <stdio.h> 21#include <math.h> 22#include <stdlib.h> 23#include <time.h> 24#ifdef _WIN32 25 26#include <process.h> 27 28#endif 29 30#include <string.h> 31#include "giants.h" 32 33 34/* definitions */ 35 36#define D 100 37#define NUM_PRIMES 6542 /* PrimePi[2^16]. */ 38 39 40/* compiler options */ 41 42#ifdef _WIN32 43#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */ 44#endif 45 46 47/* global variables */ 48 49extern giant scratch2; 50int pr[NUM_PRIMES]; 51giant xr = NULL, xs = NULL, zs = NULL, zr = NULL, xorg = NULL, 52 zorg = NULL, t1 = NULL, t2 = NULL, t3 = NULL, N = NULL, 53 gg = NULL, An = NULL, Ad = NULL; 54giant xb[D+2], zb[D+2], xzb[D+2]; 55int modmode = 0, Q, modcount = 0; 56 57 58/* function prototypes */ 59 60void ell_even(giant x1, giant z1, giant x2, giant z2, giant An, 61 giant Ad, giant N); 62void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xor, 63 giant zor, giant N); 64void ell_mul(giant xx, giant zz, int n, giant An, giant Ad, giant N); 65int least_2(int n); 66void dot(void); 67int psi_rand(); 68 69 70/************************************************************** 71 * 72 * Functions 73 * 74 **************************************************************/ 75 76 77int 78psi_rand( 79 void 80) 81{ 82 unsigned short hi; 83 unsigned short low; 84 time_t tp; 85 int result; 86 87 time(&tp); 88 low = (unsigned short)rand(); 89 hi = (unsigned short)rand(); 90 result = ((hi << 16) | low) ^ ((int)tp); 91 92 return (result & 0x7fffffff); 93} 94 95 96void 97set_random_seed( 98 void 99) 100{ 101 /* Start the random number generator at a new position. */ 102 time_t tp; 103 104 time(&tp); 105 srand((int)tp + (int)getpid()); 106} 107 108 109int 110isprime( 111 int odd 112) 113{ 114 int j; 115 int p; 116 117 for (j=1; ; j++) 118 { 119 p = pr[j]; 120 if (p*p > odd) 121 return(1); 122 if (odd % p == 0) 123 return(0); 124 } 125} 126 127 128int 129primeq( 130 int p 131) 132{ 133 register int j=3; 134 135 if ((p&1)==0) 136 return ((p==2)?1:0); 137 if (j>=p) 138 return (1); 139 while ((p%j)!=0) 140 { 141 j+=2; 142 if (j*j>p) 143 return(1); 144 } 145 return(0); 146} 147 148 149void 150s_modg( 151 giant N, 152 giant t 153) 154{ 155 ++modcount; 156 switch (modmode) 157 { 158 case 0: 159 modg(N, t); 160 break; 161 case -1: 162 mersennemod(Q, t); 163 break; 164 case 1: 165 fermatmod(Q, t); 166 break; 167 } 168} 169 170 171void 172reset_mod( 173 giant x, 174 giant N 175) 176/* Perform a divide (by the discovered factor) and switch back 177 to non-Fermat-non-Mersenne (i.e. normal) mod mode. */ 178{ 179 divg(x, N); 180 modmode = 0; 181} 182 183void 184ell_even( 185 giant x1, 186 giant z1, 187 giant x2, 188 giant z2, 189 giant An, 190 giant Ad, 191 giant N 192) 193{ 194 gtog(x1, t1); 195 addg(z1, t1); 196 squareg(t1); 197 s_modg(N, t1); 198 gtog(x1, t2); 199 subg(z1, t2); 200 squareg(t2); 201 s_modg(N, t2); 202 gtog(t1, t3); 203 subg(t2, t3); 204 gtog(t2, x2); 205 mulg(t1, x2); 206 gshiftleft(2, x2); 207 s_modg(N, x2); 208 mulg(Ad, x2); 209 s_modg(N, x2); 210 mulg(Ad, t2); 211 gshiftleft(2, t2); 212 s_modg(N, t2); 213 gtog(t3, t1); 214 mulg(An, t1); 215 s_modg(N, t1); 216 addg(t1, t2); 217 mulg(t3, t2); 218 s_modg(N, t2); 219 gtog(t2,z2); 220} 221 222 223void 224ell_odd( 225 giant x1, 226 giant z1, 227 giant x2, 228 giant z2, 229 giant xor, 230 giant zor, 231 giant N 232) 233{ 234 gtog(x1, t1); 235 subg(z1, t1); 236 gtog(x2, t2); 237 addg(z2, t2); 238 mulg(t1, t2); 239 s_modg(N, t2); 240 gtog(x1, t1); 241 addg(z1, t1); 242 gtog(x2, t3); 243 subg(z2, t3); 244 mulg(t3, t1); 245 s_modg(N, t1); 246 gtog(t2, x2); 247 addg(t1, x2); 248 squareg(x2); 249 s_modg(N, x2); 250 gtog(t2, z2); 251 subg(t1, z2); 252 squareg(z2); 253 s_modg(N, z2); 254 mulg(zor, x2); 255 s_modg(N, x2); 256 mulg(xor, z2); 257 s_modg(N, z2); 258} 259 260 261void 262ell_mul( 263 giant xx, 264 giant zz, 265 int n, 266 giant An, 267 giant Ad, 268 giant N 269) 270{ 271 unsigned int c = (unsigned int)0x80000000; 272 273 if (n==1) 274 return; 275 if (n==2) 276 { 277 ell_even(xx, zz, xx, zz, An, Ad, N); 278 return; 279 } 280 gtog(xx, xorg); 281 gtog(zz, zorg); 282 ell_even(xx, zz, xs, zs, An, Ad, N); 283 284 while((c&n) == 0) 285 { 286 c >>= 1; 287 } 288 289 c>>=1; 290 do 291 { 292 if (c&n) 293 { 294 ell_odd(xs, zs, xx, zz, xorg, zorg, N); 295 ell_even(xs, zs, xs, zs, An, Ad, N); 296 } 297 else 298 { 299 ell_odd(xx, zz, xs, zs, xorg, zorg, N); 300 ell_even(xx, zz, xx, zz, An, Ad, N); 301 } 302 c >>= 1; 303 } while(c); 304} 305 306 307 308/* From R. P. Brent, priv. comm. 1996: 309Let s > 5 be a pseudo-random seed (called $\sigma$ in the Tech. Report), 310 311 u/v = (s^2 - 5)/(4s) 312 313Then starting point is (x_1, y_1) where 314 315 x_1 = (u/v)^3 316and 317 a = (v-u)^3(3u+v)/(4u^3 v) - 2 318*/ 319 320void 321choose12( 322 giant x, 323 giant z, 324 int k, 325 giant An, 326 giant Ad, 327 giant N 328) 329{ 330 itog(k, zs); 331 gtog(zs, xs); 332 squareg(xs); 333 itog(5, t2); 334 subg(t2, xs); 335 s_modg(N, xs); 336 addg(zs, zs); 337 addg(zs, zs); 338 s_modg(N, zs); 339 gtog(xs, x); 340 squareg(x); 341 s_modg(N, x); 342 mulg(xs, x); 343 s_modg(N, x); 344 gtog(zs, z); 345 squareg(z); 346 s_modg(N, z); 347 mulg(zs, z); 348 s_modg(N, z); 349 350 /* Now for A. */ 351 gtog(zs, t2); 352 subg(xs, t2); 353 gtog(t2, t3); 354 squareg(t2); 355 s_modg(N, t2); 356 mulg(t3, t2); 357 s_modg(N, t2); /* (v-u)^3. */ 358 gtog(xs, t3); 359 addg(t3, t3); 360 addg(xs, t3); 361 addg(zs, t3); 362 s_modg(N, t3); 363 mulg(t3, t2); 364 s_modg(N, t2); /* (v-u)^3 (3u+v). */ 365 gtog(zs, t3); 366 mulg(xs, t3); 367 s_modg(N, t3); 368 squareg(xs); 369 s_modg(N, xs); 370 mulg(xs, t3); 371 s_modg(N, t3); 372 addg(t3, t3); 373 addg(t3, t3); 374 s_modg(N, t3); 375 gtog(t3, Ad); 376 gtog(t2, An); /* An/Ad is now A + 2. */ 377} 378 379 380void 381ensure( 382 int q 383) 384{ 385 int nsh, j; 386 387 N = newgiant(INFINITY); 388 if(!q) 389 { 390 gread(N,stdin); 391 q = bitlen(N) + 1; 392 } 393 nsh = q/4; /* Allowing (easily) enough space per giant, 394 since N is about 2^q, which is q bits, or 395 q/16 shorts. But squaring, etc. is allowed, 396 so we need at least q/8, and we choose q/4 397 to be conservative. */ 398 if (!xr) 399 xr = newgiant(nsh); 400 if (!zr) 401 zr = newgiant(nsh); 402 if (!xs) 403 xs = newgiant(nsh); 404 if (!zs) 405 zs = newgiant(nsh); 406 if (!xorg) 407 xorg = newgiant(nsh); 408 if (!zorg) 409 zorg = newgiant(nsh); 410 if (!t1) 411 t1 = newgiant(nsh); 412 if (!t2) 413 t2 = newgiant(nsh); 414 if (!t3) 415 t3 = newgiant(nsh); 416 if (!gg) 417 gg = newgiant(nsh); 418 if (!An) 419 An = newgiant(nsh); 420 if (!Ad) 421 Ad = newgiant(nsh); 422 for (j=0;j<D+2;j++) 423 { 424 xb[j] = newgiant(nsh); 425 zb[j] = newgiant(nsh); 426 xzb[j] = newgiant(nsh); 427 } 428} 429 430int 431bigprimeq( 432 giant x 433) 434{ 435 itog(1, t1); 436 gtog(x, t2); 437 subg(t1, t2); 438 itog(5, t1); 439 powermodg(t1, t2, x); 440 if (isone(t1)) 441 return(1); 442 return(0); 443} 444 445 446void 447dot( 448 void 449) 450{ 451 printf("."); 452 fflush(stdout); 453} 454 455/************************************************************** 456 * 457 * Main Function 458 * 459 **************************************************************/ 460 461main( 462 int argc, 463 char *argv[] 464) 465{ 466 int j, k, C, nshorts, cnt, count, 467 limitbits = 0, pass, npr, rem; 468 long B; 469 int randmode = 0; 470 471 if (!strcmp(argv[argc-1], "-r")) 472 { 473 randmode = 1; 474 if (argc > 4) 475 /* This segment only takes effect in random mode. */ 476 limitbits = atoi(argv[argc-2]); 477 } 478 else 479 { 480 randmode = 0; 481 } 482 483 modmode = 0; 484 if (argc > 2) 485 { 486 modmode = atoi(argv[1]); 487 Q = atoi(argv[2]); 488 } 489 if (modmode==0) 490 Q = 0; 491 ensure(Q); 492 if (modmode) 493 { 494 itog(1, N); 495 gshiftleft(Q, N); 496 itog(modmode, t1); 497 addg(t1, N); 498 } 499 pr[0] = 2; 500 for (k=0, npr=1;; k++) 501 { 502 if (primeq(3+2*k)) 503 { 504 pr[npr++] = 3+2*k; 505 if (npr >= NUM_PRIMES) 506 break; 507 } 508 } 509 510 if (randmode == 0) 511 { 512 printf("Sieving...\n"); 513 fflush(stdout); 514 for (j=0; j < NUM_PRIMES; j++) 515 { 516 gtog(N, t1); 517 rem = idivg(pr[j], t1); 518 if (rem == 0) 519 { 520 printf("%d ", pr[j]); 521 gtog(t1, N); 522 if (isone(N)) 523 { 524 printf("\n"); 525 exit(0); 526 } 527 else 528 { 529 printf("* "); 530 fflush(stdout); 531 } 532 --j; 533 } 534 } 535 536 if (bigprimeq(N)) 537 { 538 gout(N); 539 exit(0); 540 } 541 542 printf("\n"); 543 printf("Commencing Pollard rho...\n"); 544 fflush(stdout); 545 itog(1, gg); 546 itog(3, t1); itog(3, t2); 547 548 for (j=0; j < 15000; j++) 549 { 550 if((j%100) == 0) 551 { 552 dot(); 553 gcdg(N, gg); 554 if (!isone(gg)) 555 break; 556 } 557 squareg(t1); 558 iaddg(2, t1); 559 s_modg(N, t1); 560 squareg(t2); 561 iaddg(2, t2); 562 s_modg(N, t2); 563 squareg(t2); 564 iaddg(2, t2); 565 s_modg(N, t2); 566 gtog(t2, t3); 567 subg(t1, t3); 568 t3->sign = abs(t3->sign); 569 mulg(t3, gg); 570 s_modg(N, gg); 571 } 572 gcdg(N, gg); 573 574 if ((gcompg(N,gg) != 0) && (!isone(gg))) 575 { 576 fprintf(stdout,"\n"); 577 gout(gg); 578 reset_mod(gg, N); 579 if (isone(N)) 580 { 581 printf("\n"); 582 exit(0); 583 } 584 else 585 { 586 printf("* "); 587 fflush(stdout); 588 } 589 if (bigprimeq(N)) 590 { 591 gout(N); 592 exit(0); 593 } 594 } 595 596 printf("\n"); 597 printf("Commencing Pollard (p-1)...\n"); 598 fflush(stdout); 599 itog(1, gg); 600 itog(3, t1); 601 for (j=0; j< NUM_PRIMES; j++) 602 { 603 cnt = (int)(8*log(2.0)/log(1.0*pr[j])); 604 if (cnt < 2) 605 cnt = 1; 606 for (k=0; k< cnt; k++) 607 { 608 powermod(t1, pr[j], N); 609 } 610 itog(1, t2); 611 subg(t1, t2); 612 mulg(t2, gg); 613 s_modg(N, gg); 614 615 if (j % 100 == 0) 616 { 617 dot(); 618 gcdg(N, gg); 619 if (!isone(gg)) 620 break; 621 } 622 } 623 gcdg(N, gg); 624 if ((gcompg(N,gg) != 0) && (!isone(gg))) 625 { 626 fprintf(stdout,"\n"); 627 gout(gg); 628 reset_mod(gg, N); 629 if (isone(N)) 630 { 631 printf("\n"); 632 exit(0); 633 } 634 else 635 { 636 printf("* "); 637 fflush(stdout); 638 } 639 if (bigprimeq(N)) 640 { 641 gout(N); 642 exit(0); 643 } 644 } 645 } /* This is the end of (randmode == 0) */ 646 647 printf("\n"); 648 printf("Commencing ECM...\n"); 649 fflush(stdout); 650 651 if (randmode) 652 set_random_seed(); 653 pass = 0; 654 while (++pass) 655 { 656 if (randmode == 0) 657 { 658 if (pass <= 3) 659 { 660 B = 1000; 661 } 662 else if (pass <= 10) 663 { 664 B = 10000; 665 } 666 else if (pass <= 100) 667 { 668 B = 100000L; 669 } else 670 { 671 B = 1000000L; 672 } 673 } 674 else 675 { 676 B = 2000000L; 677 } 678 C = 50*((int)B); 679 680 /* Next, choose curve with order divisible by 16 and choose 681 * a point (xr/zr) on said curve. 682 */ 683 684 /* Order-div-12 case. 685 * cnt = 8020345; Brent's parameter for stage one discovery 686 * of 27-digit factor of F_13. 687 */ 688 689 cnt = psi_rand(); /* cnt = 8020345; */ 690 choose12(xr, zr, cnt, An, Ad, N); 691 printf("Choosing curve %d, with s = %d, B = %d, C = %d:\n", pass,cnt, B, C); fflush(stdout); 692 cnt = 0; 693 nshorts = 1; 694 count = 0; 695 for (j=0;j<nshorts;j++) 696 { 697 ell_mul(xr, zr, 1<<16, An, Ad, N); 698 ell_mul(xr, zr, 3*3*3*3*3*3*3*3*3*3*3, An, Ad, N); 699 ell_mul(xr, zr, 5*5*5*5*5*5*5, An, Ad, N); 700 ell_mul(xr, zr, 7*7*7*7*7*7, An, Ad, N); 701 ell_mul(xr, zr, 11*11*11*11, An, Ad, N); 702 ell_mul(xr, zr, 13*13*13*13, An, Ad, N); 703 ell_mul(xr, zr, 17*17*17, An, Ad, N); 704 } 705 k = 19; 706 while (k<B) 707 { 708 if (isprime(k)) 709 { 710 ell_mul(xr, zr, k, An, Ad, N); 711 if (k<100) 712 ell_mul(xr, zr, k, An, Ad, N); 713 if (cnt++ %100==0) 714 dot(); 715 } 716 k += 2; 717 } 718 count = 0; 719 720 gtog(zr, gg); 721 gcdg(N, gg); 722 if ((!isone(gg))&&(bitlen(gg)>limitbits)) 723 { 724 fprintf(stdout,"\n"); 725 gwriteln(gg, stdout); 726 fflush(stdout); 727 reset_mod(gg, N); 728 if (isone(N)) 729 { 730 printf("\n"); 731 exit(0); 732 } 733 else 734 { 735 printf("* "); 736 fflush(stdout); 737 } 738 if (bigprimeq(N)) 739 { 740 gout(N); 741 exit(0); 742 } 743 continue; 744 } 745 else 746 { 747 printf("\n"); 748 fflush(stdout); 749 } 750 751 /* Continue; Invoke, to test Stage 1 only. */ 752 k = ((int)B)/D; 753 gtog(xr, xb[0]); 754 gtog(zr, zb[0]); 755 ell_mul(xb[0], zb[0], k*D+1 , An, Ad, N); 756 gtog(xr, xb[D+1]); 757 gtog(zr, zb[D+1]); 758 ell_mul(xb[D+1], zb[D+1], (k+2)*D+1 , An, Ad, N); 759 760 for (j=1; j <= D; j++) 761 { 762 gtog(xr, xb[j]); 763 gtog(zr, zb[j]); 764 ell_mul(xb[j], zb[j], 2*j , An, Ad, N); 765 gtog(zb[j], xzb[j]); 766 mulg(xb[j], xzb[j]); 767 s_modg(N, xzb[j]); 768 } 769 modcount = 0; 770 printf("\nCommencing second stage, curve %d...\n",pass); fflush(stdout); 771 count = 0; 772 itog(1, gg); 773 774 while (1) 775 { 776 gtog(zb[0], xzb[0]); 777 mulg(xb[0], xzb[0]); 778 s_modg(N, xzb[0]); 779 mulg(zb[0], gg); 780 s_modg(N,gg); /* Accumulate. */ 781 for (j = 1; j < D; j++) 782 { 783 if (!isprime(k*D+1+ 2*j)) 784 continue; 785 786 /* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */ 787 gtog(xb[0], t1); 788 subg(xb[j], t1); 789 gtog(zb[0], t2); 790 addg(zb[j], t2); 791 mulg(t1, t2); 792 s_modg(N, t1); 793 subg(xzb[0], t2); 794 addg(xzb[j], t2); 795 s_modg(N, t2); 796 --modcount; 797 mulg(t2, gg); 798 s_modg(N, gg); 799 if((++count)%1000==0) 800 dot(); 801 } 802 803 k += 2; 804 if(k*D > C) 805 break; 806 gtog(xb[D+1], xs); 807 gtog(zb[D+1], zs); 808 ell_odd(xb[D], zb[D], xb[D+1], zb[D+1], xb[0], zb[0], N); 809 gtog(xs, xb[0]); 810 gtog(zs, zb[0]); 811 } 812 813 gcdg(N, gg); 814 if((!isone(gg))&&(bitlen(gg)>limitbits)) 815 { 816 fprintf(stdout,"\n"); 817 gwriteln(gg, stdout); 818 fflush(stdout); 819 reset_mod(gg, N); 820 if (isone(N)) 821 { 822 printf("\n"); 823 exit(0); 824 } 825 else 826 { 827 printf("* "); 828 fflush(stdout); 829 } 830 if (bigprimeq(N)) 831 { 832 gout(N); 833 exit(0); 834 } 835 continue; 836 } 837 838 printf("\n"); 839 fflush(stdout); 840 } 841 842 return 0; 843} 844 845