]> git.saurik.com Git - apple/security.git/blob - OSX/libsecurity_cryptkit/lib/CurveParamDocs/factor.c
projects / apple / security.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Security-58286.60.28.tar.gz
[apple/security.git] / OSX / libsecurity_cryptkit / lib / CurveParamDocs / factor.c
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
49 extern giant scratch2;
50 int pr[NUM_PRIMES];
51 giant 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;
54 giant xb[D+2], zb[D+2], xzb[D+2];
55 int modmode = 0, Q, modcount = 0;
56
57
58 /* function prototypes */
59
60 void ell_even(giant x1, giant z1, giant x2, giant z2, giant An,
61 giant Ad, giant N);
62 void ell_odd(giant x1, giant z1, giant x2, giant z2, giant xor,
63 giant zor, giant N);
64 void ell_mul(giant xx, giant zz, int n, giant An, giant Ad, giant N);
65 int least_2(int n);
66 void dot(void);
67 int psi_rand();
68
69
70 /**************************************************************
71 *
72 * Functions
73 *
74 **************************************************************/
75
76
77 int
78 psi_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
96 void
97 set_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
109 int
110 isprime(
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
128 int
129 primeq(
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
149 void
150 s_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
171 void
172 reset_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
183 void
184 ell_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
223 void
224 ell_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
261 void
262 ell_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:
309 Let s > 5 be a pseudo-random seed (called $\sigma$ in the Tech. Report),
310
311 u/v = (s^2 - 5)/(4s)
312
313 Then starting point is (x_1, y_1) where
314
315 x_1 = (u/v)^3
316 and
317 a = (v-u)^3(3u+v)/(4u^3 v) - 2
318 */
319
320 void
321 choose12(
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
380 void
381 ensure(
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
430 int
431 bigprimeq(
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
446 void
447 dot(
448 void
449 )
450 {
451 printf(".");
452 fflush(stdout);
453 }
454
455 /**************************************************************
456 *
457 * Main Function
458 *
459 **************************************************************/
460
461 main(
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
mirror of Security