Algorithm Alley
by Michael J. Wiener

Listing One
// To use this code, a BigInt class to handle big integers
// is needed along with the following 5 functions.

// set up an exponent to point to the top set bit for next_exponent_bit()
void init_exponent(BigInt exponent);

// get next exponent bit
unsigned long next_exponent_bit();

// return TRUE if exponent used up
int exponent_finished();

// result = result*result % modulus
void modular_square(BigInt &result, BigInt modulus);

// result = result*x % modulus
void modular_multiply(BigInt &result, BigInt x, BigInt modulus);

#define WINDOW 6
#define TABLE_LEN 32
BigInt table[TABLE_LEN];

void fill_table(BigInt message, BigInt modulus)
{
    BigInt sq;
    long i;
    sq = message;
    modular_square(sq, modulus);
    table[0] = message;
    for (i = 1; i < TABLE_LEN; i++)
    {
        table[i] = table[i - 1];
        modular_multiply(table[i], sq, modulus);
    }
}
BigInt modular_exponentiation(BigInt message, BigInt exponent, BigInt modulus)
{
    BigInt result;
    long started, num_pending_bits, num_zero_bits, i;
    unsigned long pending_bits, bit;

    fill_table(message, modulus);
    init_exponent(exponent);
    if (exponent_finished())
       return 1;
    started = 0;
    num_pending_bits = 0;
    pending_bits = 0;
    do
    {
        while (!exponent_finished() && ((bit = next_exponent_bit()) == 0))
            modular_square(result, modulus);
        num_pending_bits = pending_bits = bit;
        while (!exponent_finished() && (num_pending_bits < WINDOW))
        {
            pending_bits = (pending_bits << 1) + next_exponent_bit();
            num_pending_bits++;
        }
        if (num_pending_bits > 0)
        {
            if (!started)
            {
                if (pending_bits & 1)
                    result = table[pending_bits >> 1];
                else if (pending_bits & 2)
                {
                    result = table[pending_bits >> 2];
                    modular_square(result, modulus);
                }
                else
                {
                    result = table[(pending_bits-1) >> 1];
                    modular_multiply(result, message, modulus);
                }
                started = 1;
            }
            else
            {
                for (num_zero_bits = 0; !(pending_bits & 1); num_zero_bits++)
                    pending_bits >>= 1;
                for (i = num_zero_bits; i < num_pending_bits; i++)
                    modular_square(result, modulus);
                modular_multiply(result, table[pending_bits >> 1], modulus);
                for (i = 0; i < num_zero_bits; i++)
                    modular_square(result, modulus);
            }
        }
    } 
    while (!exponent_finished());
    return result;
}

Listing Two
// To use this code, a BigInt class to handle big integers
// is needed along with the following function.

// return true if candidate is a highly-probable prime
int miller_rabin_test(BigInt candidate);

#define SMALL_PRIME_BOUND_DIV2 2048
#define SQRT_BOUND 64
char small_prime_flags[SMALL_PRIME_BOUND_DIV2];

// sieve to find small primes
// small_prime_flag[i] == 1 means 2*i+1 is prime
void generate_small_primes()
{
    unsigned long i, sieve_val;
    small_prime_flags[0] = 0; // 1 is not prime
    for (i = 1; i < SMALL_PRIME_BOUND_DIV2; i++)
        small_prime_flags[i] = 1;
    // for each odd number, throw out its multiples
    for (sieve_val = 3; sieve_val <= SQRT_BOUND; sieve_val += 2)
        for (i = sieve_val + (sieve_val >> 1); i < SMALL_PRIME_BOUND_DIV2; 
              i += sieve_val)
            small_prime_flags[i] = 0;
}
#define SIEVE_LEN 2048
BigInt generate_large_prime(BigInt start)
{
    unsigned long small_prime, i, sp, candidate;
    char sieve_array[SIEVE_LEN];
    generate_small_primes();
    start |= 1;  // force starting point odd
    for (;; start += 2*SIEVE_LEN)
    {
        for (i = 0; i < SIEVE_LEN; i++)
            sieve_array[i] = 1;
        for (sp = 0; sp < SMALL_PRIME_BOUND_DIV2; sp++)
            if (small_prime_flags[sp])
        {
            small_prime = 2*sp + 1;  // next prime to sieve with
            // magic to find i such that small_prime divides start+2*i
            i = (small_prime - 1) - ((start - 1) % small_prime);
            if (i & 1)
                i += small_prime;
            i /= 2;
            // remove multiples of small_prime
            for (; i < SIEVE_LEN; i += small_prime)
                sieve_array[i] = 0;
        }
        // test primality of remaining candidates
        for (i = 0; i < SIEVE_LEN; i++)
            if (sieve_array[i])
        {
            candidate = start + 2*i;
            if (miller_rabin_test(candidate))
                return candidate;
       }
    }
}






3