prime.c

包含哈希,对称以及非对称的经典算法 包含经典事例

 * (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.
 *
 */
int
primeGen(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.
 */
int
primeGenStrong(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;
}