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