1/**************************************************************
2 *
3 *	fmodule.c
4 *
5 *	Factoring utilities.
6 *
7 *	Updates:
8 *		13 Apr 98   REC - creation
9 *
10 *	c. 1998 Perfectly Scientific, Inc.
11 *	All Rights Reserved.
12 *
13 *
14 *************************************************************/
15
16/* include files */
17
18#include <stdio.h>
19#include <math.h>
20#include <stdlib.h>
21#include <time.h>
22#ifdef _WIN32
23
24#include <process.h>
25
26#endif
27
28#include <string.h>
29#include "giants.h"
30#include "fmodule.h"
31#include "ellmont.h"
32
33#define NUM_PRIMES 6542 /* PrimePi[2^16]. */
34#define GENERAL_MOD 0
35#define FERMAT_MOD 1
36#define MERSENNE_MOD (-1)
37#define D 100  /* A decimation parameter for stage-2 ECM factoring. */
38
39/* compiler options */
40
41#ifdef _WIN32
42#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */
43#endif
44
45
46/* global variables */
47
48extern int pr[NUM_PRIMES];  /* Primes array from tools.c. */
49
50unsigned short factors[NUM_PRIMES], exponents[NUM_PRIMES];
51int		modmode = GENERAL_MOD;
52int		curshorts = 0;
53static giant t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL;
54static giant An = NULL, Ad = NULL;
55static point_mont pt1, pt2;
56point_mont pb[D+2];
57giant xzb[D+2];
58
59static int verbosity = 0;
60
61/**************************************************************
62 *
63 *	Functions
64 *
65 **************************************************************/
66
67int
68init_fmodule(int shorts) {
69	curshorts = shorts;
70	pb[0] = NULL;  /* To guarantee proper ECM initialization. */
71	t1 = newgiant(shorts);
72	t2 = newgiant(shorts);
73	t3 = newgiant(shorts);
74	t4 = newgiant(shorts);
75	An = newgiant(shorts);
76    Ad = newgiant(shorts);
77	pt1 = new_point_mont(shorts);
78	pt2 = new_point_mont(shorts);
79}
80
81void
82verbose(int state)
83/* Call verbose(1) for output during factoring processes,
84   call verbose(0) to silence all that.
85 */
86{
87	verbosity = state;
88}
89
90void
91dot(void)
92{
93	printf(".");
94	fflush(stdout);
95}
96
97void
98set_mod_mode(int mode)
99/* Call this with mode := 1, 0, -1, for Fermat-mod, general mod, and Mersenne mod,
100   respectively; the point being that the special cases of
101   Fermat- and Mersenne-mod are much faster than
102   general mod.  If all mods will be with respect to a number-to-be-factored,
103   of the form N = 2^m + 1, use Fermat mod; while if N = 2^m-1, use Mersenne mod.
104 */
105{
106	modmode = mode;
107}
108
109void
110special_modg(
111	giant		N,
112	giant		t
113)
114{
115	switch (modmode)
116	{
117		case MERSENNE_MOD:
118			mersennemod(bitlen(N), t);
119			break;
120		case FERMAT_MOD:
121			fermatmod(bitlen(N)-1, t);
122			break;
123	    default:
124			modg(N, t);
125			break;
126	}
127}
128
129unsigned short *
130prime_list() {
131	return(&factors[0]);
132}
133
134unsigned short *
135exponent_list() {
136	return(&exponents[0]);
137}
138
139int
140sieve(giant N, int sievelim)
141/* Returns number of N's prime factors < min(sievelim, 2^16),
142   with N reduced accordingly by said factors.
143   The n-th entry of factors[] becomes the n-th prime
144   factor of N, with corresponding exponent
145   becoming the n-th element of exponents[].
146 */
147{	int j, pcount, rem;
148	unsigned short pri;
149
150		pcount = 0;
151		exponents[0] = 0;
152		for (j=0; j < NUM_PRIMES; j++)
153		{
154			pri = pr[j];
155			if(pri > sievelim) break;
156			do {
157				gtog(N, t1);
158				rem = idivg(pri, t1);
159				if(rem == 0) {
160					++exponents[pcount];
161					gtog(t1, N);
162				}
163			} while(rem == 0);
164			if(exponents[pcount] > 0) {
165				factors[pcount] = pr[j];
166				++pcount;
167				exponents[pcount] = 0;
168			}
169		}
170		return(pcount);
171}
172
173int
174pollard_rho(giant N, giant fact, int steps, int abort)
175/* Perform Pollard-rho x:= 3; loop(x:= x^2 + 2), a total of steps times.
176   Parameter fact will be a nontrivial factor found, in which case
177   N is also modified as: N := N/fact.
178   The function returns 0 if no nontrivial factor found, else returns 1.
179   The abort parameter, when set, causes the factorer to exit on the
180   first nontrivial factor found (the requisite GCD is checked
181   every 1000 steps).  If abort := 0, the full number
182   of steps are always performed, then one solitary GCD is taken,
183   before exit.
184 */
185{
186	int j, found = 0;
187
188	itog(3, t1);
189	gtog(t1, t2);
190	itog(1, fact);
191	for(j=0; j < steps; j++) {
192		squareg(t1); iaddg(2, t1); special_modg(N, t1);
193		squareg(t2); iaddg(2, t2); special_modg(N, t2);
194		squareg(t2); iaddg(2, t2); special_modg(N, t2);
195		gtog(t2, t3); subg(t1,t3); mulg(t3, fact); special_modg(N, fact);
196		if(((j % 1000 == 999) && abort) || (j == steps-1)) {
197			if(verbosity) dot();
198			gcdg(N, fact);
199			if(!isone(fact)) {
200				found = (gcompg(N, fact) == 0) ? 0 : 1;
201				break;
202			}
203		}
204	}
205    if(verbosity) { printf("\n"); fflush(stdout); }
206	if(found) {
207		divg(fact, N);
208		return(1);
209	}
210	itog(1, fact);
211	return(0);
212}
213
214int
215pollard_pminus1(giant N, giant fact, int lim, int abort)
216/* Perform Pollard-(p-1); where we test
217
218		GCD[N, 3^P - 1],
219
220   where P is an accumulation of primes <= min(lim, 2^16),
221   to appropriate powers.
222   Parameter fact will be a nontrivial factor found, in which case
223   N is also modified as: N := N/fact.
224   The function returns 0 if no nontrivial factor found, else returns 1.
225   The abort parameter, when set, causes the factorer to exit on the
226   first nontrivial factor found (the requisite GCD is checked
227   every 100 steps).  If abort := 0, the full number
228   of steps are always performed, then one solitary GCD is taken,
229   before exit.
230 */
231{  int cnt, j, k, pri, found = 0;
232
233   itog(3, fact);
234   for (j=0; j< NUM_PRIMES; j++)
235	{
236			pri = pr[j];
237			if((pri > lim) || (j == NUM_PRIMES-1) || (abort && (j % 100 == 99))) {
238				if(verbosity) dot();
239				gtog(fact, t1);
240				itog(1, t2);
241				subg(t2, t1);
242				special_modg(N, t1);
243				gcdg(N, t1);
244				if(!isone(t1)) {
245					found = (gcompg(N, t1) == 0) ? 0 : 1;
246					break;
247				}
248				if(pri > lim) break;
249            }
250			if(pri < 19) { cnt = 20-pri;  /* Smaller primes get higher powers. */
251			} else if(pri < 100) {
252						cnt = 2;
253				   } else cnt = 1;
254			for (k=0; k< cnt; k++)
255			{
256				powermod(fact, pri, N);
257			}
258	}
259    if(verbosity) { printf("\n"); fflush(stdout); }
260	if(found) {
261		gtog(t1, fact);
262		divg(fact, N);
263		return(1);
264	}
265	itog(1, fact);
266	return(0);
267}
268
269int
270ecm(giant N, giant fact, int S, unsigned int B, unsigned int C)
271/* Perform elliptic curve method (ECM), with:
272   Brent seed parameter = S
273   Stage-one limit = B
274   Stage-two limit = C
275   This function:
276		returns 1 if nontrivial factor is found in stage 1 of ECM;
277		returns 2 if nontrivial factor is found in stage 2 of ECM;
278		returns 0 otherwise.
279   In the positive return, parameter fact is the factor and N := N/fact.
280 */
281{	unsigned int pri, q;
282	int j, cnt, count, k;
283
284    if(verbosity) {
285	  	printf("Finding curve and point, B = %u, C = %u, seed = %d...", B, C, S);
286	  	fflush(stdout);
287    }
288	find_curve_point_brent(pt1, S, An, Ad, N);
289    if(verbosity) {
290	  	printf("\n");
291		printf("Commencing stage 1 of ECM...\n");
292		fflush(stdout);
293    }
294
295	q = pr[NUM_PRIMES-1];
296	count = 0;
297	for(j=0; ; j++) {
298		if(j < NUM_PRIMES) {
299			pri = pr[j];
300		} else {
301			q += 2;
302			if(primeq(q)) pri = q;
303				else continue;
304		}
305		if(verbosity) if((++count) % 100 == 0) dot();
306		if(pri > B) break;
307		if(pri < 19) { cnt = 20-pri;
308			} else if(pri < 100) {
309						cnt = 2;
310				   } else cnt = 1;
311		for(k = 0; k < cnt; k++)
312			ell_mul_int_brent(pt1, pri, An, Ad, N);
313	}
314    k = 19;
315	while (k<B)
316	{
317		if (primeq(k))
318		{
319			ell_mul_int_brent(pt1, k, An, Ad, N);
320			if (k<100)
321					ell_mul_int_brent(pt1, k, An, Ad, N);
322			if (cnt++ %100==0)
323					dot();
324		}
325		k += 2;
326	}
327    if(verbosity) { printf("\n"); fflush(stdout); }
328
329/* Next, test stage-1 attempt. */
330    gtog(pt1->z, fact);
331	gcdg(N, fact);
332    if((!isone(fact)) && (gcompg(N, fact) != 0)) {
333		divg(fact, N);
334		return(1);
335	}
336	if(B >= C) {  /* No stage 2 planned. */
337		itog(1, fact);
338		return(0);
339	}
340
341/* Commence second stage of ECM. */
342    if(verbosity) {
343	  	printf("\n");
344		printf("Commencing stage 2 of ECM...\n");
345		fflush(stdout);
346    }
347	if(pb[0] == NULL) {
348		for(k=0; k < D+2; k++) {
349				pb[k] = new_point_mont(curshorts);
350				xzb[k] = newgiant(curshorts);
351
352		}
353	}
354	k = ((int)B)/D;
355	ptop_mont(pt1, pb[0]);
356	ell_mul_int_brent(pb[0], k*D+1 , An, Ad, N);
357	ptop_mont(pt1, pb[D+1]);
358	ell_mul_int_brent(pb[D+1], (k+2)*D+1 , An, Ad, N);
359
360	for (j=1; j <= D; j++)
361	{
362		ptop_mont(pt1, pb[j]);
363		ell_mul_int_brent(pb[j], 2*j , An, Ad, N);
364		gtog(pb[j]->z, xzb[j]);
365		mulg(pb[j]->x, xzb[j]);
366		special_modg(N, xzb[j]);
367	}
368	itog(1, fact);
369	count = 0;
370	while (1) {
371		if(verbosity) if((++count) % 10 == 0) dot();
372		gtog(pb[0]->z, xzb[0]);
373		mulg(pb[0]->x, xzb[0]);
374		special_modg(N, xzb[0]);
375		mulg(pb[0]->z, fact);
376		special_modg(N, fact); /* Accumulate. */
377		for (j = 1; j < D; j++) {
378				if (!primeq(k*D+1+2*j)) continue;
379/* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */
380				gtog(pb[0]->x, t1);
381				subg(pb[j]->x, t1);
382				gtog(pb[0]->z, t2);
383				addg(pb[j]->z, t2);
384				mulg(t1, t2);
385				special_modg(N, t1);
386				subg(xzb[0], t2);
387				addg(xzb[j], t2);
388				special_modg(N, t2);
389				mulg(t2, fact);
390				special_modg(N, fact);
391		}
392		k += 2;
393		if(k*D > C)
394		break;
395		ptop_mont(pb[D+1], pt2);
396		ell_odd_brent(pb[D], pb[D+1], pb[0], N);
397		ptop_mont(pt2, pb[0]);
398	}
399    if(verbosity) { printf("\n"); fflush(stdout); }
400
401	gcdg(N, fact);
402    if((!isone(fact)) && (gcompg(N, fact) != 0)) {
403		divg(fact, N);
404		return(2);  /* Return value of 2 for stage-2 success! */
405	}
406	itog(1, fact);
407	return(0);
408}
409
410
411