[apple/security.git] / OSX / include / security_cryptkit / CurveParamDocs / fmodule.c

/**************************************************************
 *
 *	fmodule.c
 *
 *	Factoring utilities.
 *
 *	Updates:
 *		13 Apr 98   REC - creation
 *
 *	c. 1998 Perfectly Scientific, Inc.
 *	All Rights Reserved.
 *
 *
 *************************************************************/

/* include files */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#ifdef _WIN32 

#include <process.h>

#endif

#include <string.h>
#include "giants.h"
#include "fmodule.h"
#include "ellmont.h"

#define NUM_PRIMES 6542 /* PrimePi[2^16]. */
#define GENERAL_MOD 0
#define FERMAT_MOD 1
#define MERSENNE_MOD (-1)
#define D 100  /* A decimation parameter for stage-2 ECM factoring. */

/* compiler options */

#ifdef _WIN32
#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */
#endif


/* global variables */

extern int pr[NUM_PRIMES];  /* Primes array from tools.c. */

unsigned short factors[NUM_PRIMES], exponents[NUM_PRIMES];
int		modmode = GENERAL_MOD;
int		curshorts = 0;
static giant t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL;
static giant An = NULL, Ad = NULL;
static point_mont pt1, pt2;
point_mont pb[D+2];
giant xzb[D+2];

static int verbosity = 0;

/**************************************************************
 *
 *	Functions
 *
 **************************************************************/

int
init_fmodule(int shorts) {
	curshorts = shorts;
	pb[0] = NULL;  /* To guarantee proper ECM initialization. */
	t1 = newgiant(shorts);
	t2 = newgiant(shorts);
	t3 = newgiant(shorts);
	t4 = newgiant(shorts);
	An = newgiant(shorts);
    Ad = newgiant(shorts);
	pt1 = new_point_mont(shorts);
	pt2 = new_point_mont(shorts);
}

void
verbose(int state)
/* Call verbose(1) for output during factoring processes, 
   call verbose(0) to silence all that.
 */ 
{
	verbosity = state;
}

void
dot(void)
{
	printf(".");
	fflush(stdout);
}

void
set_mod_mode(int mode) 
/* Call this with mode := 1, 0, -1, for Fermat-mod, general mod, and Mersenne mod, 
   respectively; the point being that the special cases of
   Fermat- and Mersenne-mod are much faster than
   general mod.  If all mods will be with respect to a number-to-be-factored,
   of the form N = 2^m + 1, use Fermat mod; while if N = 2^m-1, use Mersenne mod.
 */
{
	modmode = mode;
}

void
special_modg(
	giant		N,
	giant		t
)
{
	switch (modmode)
	{
		case MERSENNE_MOD:
			mersennemod(bitlen(N), t);
			break;
		case FERMAT_MOD:
			fermatmod(bitlen(N)-1, t);
			break;
	    default:
			modg(N, t);
			break;
	}
}

unsigned short *
prime_list() {
	return(&factors[0]);
}

unsigned short *
exponent_list() {
	return(&exponents[0]);
}

int
sieve(giant N, int sievelim)
/* Returns number of N's prime factors < min(sievelim, 2^16),
   with N reduced accordingly by said factors.
   The n-th entry of factors[] becomes the n-th prime
   factor of N, with corresponding exponent
   becoming the n-th element of exponents[]. 
 */
{	int j, pcount, rem;
	unsigned short pri;

		pcount = 0;
		exponents[0] = 0;
		for (j=0; j < NUM_PRIMES; j++)
		{
			pri = pr[j];
			if(pri > sievelim) break;
			do {
				gtog(N, t1);
				rem = idivg(pri, t1);
				if(rem == 0) {
					++exponents[pcount];
					gtog(t1, N);
				}
			} while(rem == 0);
			if(exponents[pcount] > 0) {
				factors[pcount] = pr[j];
				++pcount;
				exponents[pcount] = 0;
			}
		}
		return(pcount);
}

int
pollard_rho(giant N, giant fact, int steps, int abort)
/* Perform Pollard-rho x:= 3; loop(x:= x^2 + 2), a total of steps times.
   Parameter fact will be a nontrivial factor found, in which case
   N is also modified as: N := N/fact.
   The function returns 0 if no nontrivial factor found, else returns 1.
   The abort parameter, when set, causes the factorer to exit on the
   first nontrivial factor found (the requisite GCD is checked
   every 1000 steps).  If abort := 0, the full number
   of steps are always performed, then one solitary GCD is taken,
   before exit.
 */
{	
	int j, found = 0;

	itog(3, t1);
	gtog(t1, t2);
	itog(1, fact);
	for(j=0; j < steps; j++) {
		squareg(t1); iaddg(2, t1); special_modg(N, t1);
		squareg(t2); iaddg(2, t2); special_modg(N, t2);
		squareg(t2); iaddg(2, t2); special_modg(N, t2);
		gtog(t2, t3); subg(t1,t3); mulg(t3, fact); special_modg(N, fact);
		if(((j % 1000 == 999) && abort) || (j == steps-1)) {
			if(verbosity) dot();
			gcdg(N, fact);
			if(!isone(fact)) {
				found = (gcompg(N, fact) == 0) ? 0 : 1;
				break;
			}
		}			
	}
    if(verbosity) { printf("\n"); fflush(stdout); }
	if(found) {
		divg(fact, N);
		return(1);
	}
	itog(1, fact);
	return(0);		
}

int
pollard_pminus1(giant N, giant fact, int lim, int abort)
/* Perform Pollard-(p-1); where we test
	
		GCD[N, 3^P - 1],

   where P is an accumulation of primes <= min(lim, 2^16), 
   to appropriate powers.
   Parameter fact will be a nontrivial factor found, in which case
   N is also modified as: N := N/fact.
   The function returns 0 if no nontrivial factor found, else returns 1.
   The abort parameter, when set, causes the factorer to exit on the
   first nontrivial factor found (the requisite GCD is checked
   every 100 steps).  If abort := 0, the full number
   of steps are always performed, then one solitary GCD is taken,
   before exit.  
 */
{  int cnt, j, k, pri, found = 0;

   itog(3, fact);
   for (j=0; j< NUM_PRIMES; j++)
	{
			pri = pr[j];
			if((pri > lim) || (j == NUM_PRIMES-1) || (abort && (j % 100 == 99))) {
				if(verbosity) dot();
				gtog(fact, t1);
				itog(1, t2);
				subg(t2, t1);
				special_modg(N, t1);
				gcdg(N, t1);
				if(!isone(t1)) {
					found = (gcompg(N, t1) == 0) ? 0 : 1;
					break;
				}
				if(pri > lim) break;
            }
			if(pri < 19) { cnt = 20-pri;  /* Smaller primes get higher powers. */
			} else if(pri < 100) {
						cnt = 2; 
				   } else cnt = 1;
			for (k=0; k< cnt; k++)
			{
				powermod(fact, pri, N);
			}
	}
    if(verbosity) { printf("\n"); fflush(stdout); }
	if(found) {
		gtog(t1, fact);
		divg(fact, N);
		return(1);
	}
	itog(1, fact);
	return(0);		
}

int
ecm(giant N, giant fact, int S, unsigned int B, unsigned int C)
/* Perform elliptic curve method (ECM), with:
   Brent seed parameter = S
   Stage-one limit = B
   Stage-two limit = C
   This function:
		returns 1 if nontrivial factor is found in stage 1 of ECM;
		returns 2 if nontrivial factor is found in stage 2 of ECM;
		returns 0 otherwise.
   In the positive return, parameter fact is the factor and N := N/fact.
 */
{	unsigned int pri, q;
	int j, cnt, count, k;

    if(verbosity) {
	  	printf("Finding curve and point, B = %u, C = %u, seed = %d...", B, C, S);
	  	fflush(stdout);
    }
	find_curve_point_brent(pt1, S, An, Ad, N);
    if(verbosity) {
	  	printf("\n"); 
		printf("Commencing stage 1 of ECM...\n");
		fflush(stdout);
    }

	q = pr[NUM_PRIMES-1];
	count = 0;
	for(j=0; ; j++) {
		if(j < NUM_PRIMES) {
			pri = pr[j];
		} else {
			q += 2;
			if(primeq(q)) pri = q; 
				else continue;
		}
		if(verbosity) if((++count) % 100 == 0) dot();
		if(pri > B) break;
		if(pri < 19) { cnt = 20-pri; 
			} else if(pri < 100) {
						cnt = 2; 
				   } else cnt = 1;
		for(k = 0; k < cnt; k++)	
			ell_mul_int_brent(pt1, pri, An, Ad, N); 
	}
    k = 19;
	while (k<B)
	{
		if (primeq(k))
		{
			ell_mul_int_brent(pt1, k, An, Ad, N);
			if (k<100)
					ell_mul_int_brent(pt1, k, An, Ad, N);
			if (cnt++ %100==0)
					dot();
		}
		k += 2;
	}
    if(verbosity) { printf("\n"); fflush(stdout); }

/* Next, test stage-1 attempt. */	
    gtog(pt1->z, fact);
	gcdg(N, fact);
    if((!isone(fact)) && (gcompg(N, fact) != 0)) {
		divg(fact, N);
		return(1);
	}
	if(B >= C) {  /* No stage 2 planned. */
		itog(1, fact);
		return(0);
	}

/* Commence second stage of ECM. */
    if(verbosity) {
	  	printf("\n"); 
		printf("Commencing stage 2 of ECM...\n");
		fflush(stdout);
    }
	if(pb[0] == NULL) {
		for(k=0; k < D+2; k++) {
				pb[k] = new_point_mont(curshorts);
				xzb[k] = newgiant(curshorts);

		}
	}
	k = ((int)B)/D;
	ptop_mont(pt1, pb[0]);
	ell_mul_int_brent(pb[0], k*D+1 , An, Ad, N);
	ptop_mont(pt1, pb[D+1]);
	ell_mul_int_brent(pb[D+1], (k+2)*D+1 , An, Ad, N);

	for (j=1; j <= D; j++)
	{
		ptop_mont(pt1, pb[j]);
		ell_mul_int_brent(pb[j], 2*j , An, Ad, N);
		gtog(pb[j]->z, xzb[j]);
		mulg(pb[j]->x, xzb[j]);
		special_modg(N, xzb[j]);
	}
	itog(1, fact);
	count = 0;
	while (1) {
		if(verbosity) if((++count) % 10 == 0) dot();
		gtog(pb[0]->z, xzb[0]);
		mulg(pb[0]->x, xzb[0]);
		special_modg(N, xzb[0]);
		mulg(pb[0]->z, fact);
		special_modg(N, fact); /* Accumulate. */
		for (j = 1; j < D; j++) {
				if (!primeq(k*D+1+2*j)) continue;
/* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */
				gtog(pb[0]->x, t1);
				subg(pb[j]->x, t1);
				gtog(pb[0]->z, t2);
				addg(pb[j]->z, t2);
				mulg(t1, t2);
				special_modg(N, t1);
				subg(xzb[0], t2);
				addg(xzb[j], t2);
				special_modg(N, t2);
				mulg(t2, fact);
				special_modg(N, fact);
		}
		k += 2;
		if(k*D > C)
		break;
		ptop_mont(pb[D+1], pt2);
		ell_odd_brent(pb[D], pb[D+1], pb[0], N);
		ptop_mont(pt2, pb[0]);
	}
    if(verbosity) { printf("\n"); fflush(stdout); }

	gcdg(N, fact);
    if((!isone(fact)) && (gcompg(N, fact) != 0)) {
		divg(fact, N);
		return(2);  /* Return value of 2 for stage-2 success! */
	}
	itog(1, fact);
	return(0);	
}
Commit	Line	Data
b1ab9ed8 A	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
	48	extern int pr[NUM_PRIMES]; /* Primes array from tools.c. */
	49
	50	unsigned short factors[NUM_PRIMES], exponents[NUM_PRIMES];
	51	int modmode = GENERAL_MOD;
	52	int curshorts = 0;
	53	static giant t1 = NULL, t2 = NULL, t3 = NULL, t4 = NULL;
	54	static giant An = NULL, Ad = NULL;
	55	static point_mont pt1, pt2;
	56	point_mont pb[D+2];
	57	giant xzb[D+2];
	58
	59	static int verbosity = 0;
	60
	61	/**************************************************************
	62	*
	63	* Functions
	64	*
65	**************************************************************/
66
67	int
68	init_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
81	void
82	verbose(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
90	void
91	dot(void)
92	{
93	printf(".");
94	fflush(stdout);
95	}
96
97	void
98	set_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
109	void
110	special_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
129	unsigned short *
130	prime_list() {
131	return(&factors[0]);
132	}
133
134	unsigned short *
135	exponent_list() {
136	return(&exponents[0]);
137	}
138
139	int
140	sieve(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
173	int
174	pollard_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
214	int
215	pollard_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
269	int
270	ecm(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(kD+1+2j)) 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