* (x^e)^d == x (mod p), then d*e == 1 (mod p-1), so gcd(e,p-1) must be 1. * The prime returned is constrained to not be congruent to 1 modulo * any of the zero-terminated list of 16-bit numbers.  Note that this list * should contain all the small prime factors of e.  (You'll have to test * for large prime factors of e elsewhere, but the chances of needing to * generate another prime are low.) * * The list is terminated by a 0, and may be empty. * */intprimeGen(struct BigNum *bn, unsigned (*rand)(unsigned),         int (*f)(void *arg, int c), void *arg, unsigned exponent, ...){	int retval;	int modexps = 0;	unsigned short offsets[SHUFFLE];	unsigned i, j;	unsigned p, q, prev;	struct BigNum a, e;#ifdef MSDOS	unsigned char *sieve;#else	unsigned char sieve[SIEVE];#endif#ifdef MSDOS	sieve = lbnMemAlloc(SIEVE);	if (!sieve)		return -1;#endif	bnBegin(&a);	bnBegin(&e);#if 0	/* Self-test (not used for production) */{	struct BigNum t;	static unsigned char const prime1[] = {5};	static unsigned char const prime2[] = {7};	static unsigned char const prime3[] = {11};	static unsigned char const prime4[] = {1, 1}; /* 257 */	static unsigned char const prime5[] = {0xFF, 0xF1}; /* 65521 */	static unsigned char const prime6[] = {1, 0, 1}; /* 65537 */	static unsigned char const prime7[] = {1, 0, 3}; /* 65539 */	/* A small prime: 1234567891 */	static unsigned char const prime8[] = {0x49, 0x96, 0x02, 0xD3};	/* A slightly larger prime: 12345678901234567891 */	static unsigned char const prime9[] = {		0xAB, 0x54, 0xA9, 0x8C, 0xEB, 0x1F, 0x0A, 0xD3 };	/*	 * No, 123456789012345678901234567891 isn't prime; it's just a	 * lucky, easy-to-remember conicidence.  (You have to go to	 * ...4567907 for a prime.)	 */	static struct {		unsigned char const *prime;		unsigned size;	} const primelist[] = {		{ prime1, sizeof(prime1) },		{ prime2, sizeof(prime2) },		{ prime3, sizeof(prime3) },		{ prime4, sizeof(prime4) },		{ prime5, sizeof(prime5) },		{ prime6, sizeof(prime6) },		{ prime7, sizeof(prime7) },		{ prime8, sizeof(prime8) },		{ prime9, sizeof(prime9) } };	bnBegin(&t);	for (i = 0; i < sizeof(primelist)/sizeof(primelist[0]); i++) {			bnInsertBytes(&t, primelist[i].prime, 0,				      primelist[i].size);			bnCopy(&e, &t);			(void)bnSubQ(&e, 1);			bnTwoExpMod(&a, &e, &t);			p = bnBits(&a);			if (p != 1) {				printf(			"Bug: Fermat(2) %u-bit output (1 expected)\n", p);				fputs("Prime = 0x", stdout);				for (j = 0; j < primelist[i].size; j++)					printf("%02X", primelist[i].prime[j]);				putchar('\n');			}			bnSetQ(&a, 3);			bnExpMod(&a, &a, &e, &t);			if (p != 1) {				printf(			"Bug: Fermat(3) %u-bit output (1 expected)\n", p);				fputs("Prime = 0x", stdout);				for (j = 0; j < primelist[i].size; j++)					printf("%02X", primelist[i].prime[j]);				putchar('\n');			}		}	bnEnd(&t);}#endif	/* First, make sure that bn is odd. */	if ((bnLSWord(bn) & 1) == 0)		(void)bnAddQ(bn, 1);retry:	/* Then build a sieve starting at bn. */	sieveBuild(sieve, SIEVE, bn, 2, 0);	/* Do the extra exponent sieving */	if (exponent) {		va_list ap;		unsigned t = exponent;		va_start(ap, exponent);		do {			/* The exponent had better be odd! */			assert(t & 1);			i = bnModQ(bn, t);			/* Find 1-i */			if (i == 0)				i = 1;			else if (--i)				i = t - i;			/* Divide by 2, modulo the exponent */			i = (i & 1) ? i/2 + t/2 + 1 : i/2;			/* Remove all following multiples from the sieve. */			sieveSingle(sieve, SIEVE, i, t);			/* Get the next exponent value */			t = va_arg(ap, unsigned);		} while (t);		va_end(ap);	}	/* Fill up the offsets array with the first SHUFFLE candidates */	i = p = 0;	/* Get first prime */	if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0)		offsets[i++] = p;	/*	 * If we have a prime and we have a shuffle function so it makes	 * to fill up the table, so do.	 */	if (rand && i) {		do {			p = sieveSearch(sieve, SIEVE, p);			if (p == 0)				break;			offsets[i++] = p;		} while (i < SHUFFLE);	}	/* Choose a random candidate for experimentation */	prev = 0;	while (i) {		/* Pick a random entry from the shuffle table */		j = rand ? rand(i) : 0;		q = offsets[j];		/* Replace the entry with some more data, if possible */		if (p && (p = sieveSearch(sieve, SIEVE, p)) != 0)			offsets[j] = p;		else			offsets[j] = offsets[--i];		/* Adjust bn to have the right value */		if (q > prev) {			if (bnAddQ(bn, q-prev) < 0)				goto failed;			if (bnAddQ(bn, q-prev) < 0)				goto failed;		} else {			if (bnSubQ(bn, prev-q))				goto failed;			if (bnSubQ(bn, prev-q))				goto failed;		}		prev = q;		/* Now do the Fermat tests */		retval = primeTest(bn, &e, &a, f, arg);		if (retval <= 0)			goto done;	/* Success or error */		modexps += retval;		if (f && (retval = f(arg, '.')) < 0)			goto done;	}	/* Ran out of sieve space - increase bn and keep trying. */	if (bnSubQ(bn, (SIEVE*8)-prev))		goto failed;	if (bnSubQ(bn, (SIEVE*8)-prev))		goto failed;	if (f && (retval = f(arg, '/')) < 0)		goto done;	goto retry;failed:	retval = -1;done:	bnEnd(&e);	bnEnd(&a);	lbnMemWipe(offsets, sizeof(offsets));#ifdef MSDOS	lbnMemFree(sieve, SIEVE);#else	lbnMemWipe(sieve, sizeof(sieve));#endif	return retval < 0 ? retval : modexps;}/* * Similar, but searches forward from the given starting value in steps of * "step" rather than 1.  The step size must be even, and bn must be odd. * Among other possibilities, this can be used to generate "strong" * primes, where p-1 has a large prime factor. */intprimeGenStrong(struct BigNum *bn, struct BigNum const *step,	int (*f)(void *arg, int c), void *arg){	int retval;	unsigned p, prev;	struct BigNum a, e;	int modexps = 0;#ifdef MSDOS	unsigned char *sieve;#else	unsigned char sieve[SIEVE];#endif#ifdef MSDOS	sieve = lbnMemAlloc(SIEVE);	if (!sieve)		return -1;#endif	/* Step must be even and bn must be odd */	assert((bnLSWord(step) & 1) == 0);	assert((bnLSWord(bn) & 1) == 1);	bnBegin(&a);	bnBegin(&e);	for (;;) {		if (sieveBuildBig(sieve, SIEVE, bn, step, 0) < 0)			goto failed;		p = prev = 0;		if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) {			do {				/*				 * Adjust bn to have the right value,				 * adding (p-prev) * 2*step.				 */				assert(p >= prev);				/* Compute delta into a */				if (bnMulQ(&a, step, p-prev) < 0)					goto failed;				if (bnAdd(bn, &a) < 0)					goto failed;				prev = p;				retval = primeTest(bn, &e, &a, f, arg);				if (retval <= 0)					goto done;	/* Success! */				modexps += retval;				if (f && (retval = f(arg, '.')) < 0)					goto done;				/* And try again */				p = sieveSearch(sieve, SIEVE, p);			} while (p);		}		/* Ran out of sieve space - increase bn and keep trying. */#if SIEVE*8 == 65536		/* Corner case that will never actually happen */		if (!prev) {			if (bnAdd(bn, step) < 0)				goto failed;			p = 65535;		} else {			p = (unsigned)(SIEVE*8 - prev);		}#else		p = SIEVE*8 - prev;#endif		if (bnMulQ(&a, step, p) < 0 || bnAdd(bn, &a) < 0)			goto failed;		if (f && (retval = f(arg, '/')) < 0)			goto done;	} /* for (;;) */failed:	retval = -1;done:	bnEnd(&e);	bnEnd(&a);#ifdef MSDOS	lbnMemFree(sieve, SIEVE);#else	lbnMemWipe(sieve, sizeof(sieve));#endif	return retval < 0 ? retval : modexps;}