if (factor) multipliers[NPRIMES] = bignum_mod_short(factor, modulus); else multipliers[NPRIMES] = 1;
Bignum primegen(int bits, int modulus, int residue, Bignum factor,		int phase, progfn_t pfn, void *pfnparam){    int i, k, v, byte, bitsleft, check, checks;    unsigned long delta;    unsigned long moduli[NPRIMES + 1];    unsigned long residues[NPRIMES + 1];    unsigned long multipliers[NPRIMES + 1];    Bignum p, pm1, q, wqp, wqp2;    int progress = 0;    byte = 0;    bitsleft = 0;  STARTOVER:    pfn(pfnparam, PROGFN_PROGRESS, phase, ++progress);    /*     * Generate a k-bit random number with top and bottom bits set.     * Alternatively, if `factor' is nonzero, generate a k-bit     * random number with the top bit set and the bottom bit clear,     * multiply it by `factor', and add one.     */    p = bn_power_2(bits - 1);    for (i = 0; i < bits; i++) {	if (i == 0 || i == bits - 1)	    v = (i != 0 || !factor) ? 1 : 0;	else {	    if (bitsleft <= 0)		bitsleft = 8, byte = random_byte();	    v = byte & 1;	    byte >>= 1;	    bitsleft--;	}	bignum_set_bit(p, i, v);    }    if (factor) {	Bignum tmp = p;	p = bigmul(tmp, factor);	freebn(tmp);	assert(bignum_bit(p, 0) == 0);	bignum_set_bit(p, 0, 1);    }    /*     * Ensure this random number is coprime to the first few     * primes, by repeatedly adding either 2 or 2*factor to it     * until it is.     */    for (i = 0; i < NPRIMES; i++) {	moduli[i] = primes[i];	residues[i] = bignum_mod_short(p, primes[i]);	if (factor)	    multipliers[i] = bignum_mod_short(factor, primes[i]);	else	    multipliers[i] = 1;    }    moduli[NPRIMES] = modulus;    residues[NPRIMES] = (bignum_mod_short(p, (unsigned short) modulus)			 + modulus - residue);    delta = 0;    while (1) {	for (i = 0; i < (sizeof(moduli) / sizeof(*moduli)); i++)	    if (!((residues[i] + delta * multipliers[i]) % moduli[i]))		break;	if (i < (sizeof(moduli) / sizeof(*moduli))) {	/* we broke */	    delta += 2;	    if (delta > 65536) {		freebn(p);		goto STARTOVER;	    }	    continue;	}	break;    }    q = p;    if (factor) {	Bignum tmp;	tmp = bignum_from_long(delta);	p = bigmuladd(tmp, factor, q);	freebn(tmp);    } else {	p = bignum_add_long(q, delta);    }    freebn(q);    /*     * Now apply the Miller-Rabin primality test a few times. First     * work out how many checks are needed.     */    checks = 27;    if (bits >= 150)	checks = 18;    if (bits >= 200)	checks = 15;    if (bits >= 250)	checks = 12;    if (bits >= 300)	checks = 9;    if (bits >= 350)	checks = 8;    if (bits >= 400)	checks = 7;    if (bits >= 450)	checks = 6;    if (bits >= 550)	checks = 5;    if (bits >= 650)	checks = 4;    if (bits >= 850)	checks = 3;    if (bits >= 1300)	checks = 2;    /*     * Next, write p-1 as q*2^k.     */    for (k = 0; bignum_bit(p, k) == !k; k++)	continue;	/* find first 1 bit in p-1 */    q = bignum_rshift(p, k);    /* And store p-1 itself, which we'll need. */    pm1 = copybn(p);    decbn(pm1);    /*     * Now, for each check ...     */    for (check = 0; check < checks; check++) {	Bignum w;	/*	 * Invent a random number between 1 and p-1 inclusive.	 */	while (1) {	    w = bn_power_2(bits - 1);	    for (i = 0; i < bits; i++) {		if (bitsleft <= 0)		    bitsleft = 8, byte = random_byte();		v = byte & 1;		byte >>= 1;		bitsleft--;		bignum_set_bit(w, i, v);	    }	    bn_restore_invariant(w);	    if (bignum_cmp(w, p) >= 0 || bignum_cmp(w, Zero) == 0) {		freebn(w);		continue;	    }	    break;	}	pfn(pfnparam, PROGFN_PROGRESS, phase, ++progress);	/*	 * Compute w^q mod p.	 */	wqp = modpow(w, q, p);	freebn(w);	/*	 * See if this is 1, or if it is -1, or if it becomes -1	 * when squared at most k-1 times.	 */	if (bignum_cmp(wqp, One) == 0 || bignum_cmp(wqp, pm1) == 0) {	    freebn(wqp);	    continue;	}	for (i = 0; i < k - 1; i++) {	    wqp2 = modmul(wqp, wqp, p);	    freebn(wqp);	    wqp = wqp2;	    if (bignum_cmp(wqp, pm1) == 0)		break;	}	if (i < k - 1) {	    freebn(wqp);	    continue;	}	/*	 * It didn't. Therefore, w is a witness for the	 * compositeness of p.	 */	freebn(p);	freebn(pm1);	freebn(q);	goto STARTOVER;    }    /*     * We have a prime!     */    freebn(q);    freebn(pm1);    return p;}