How can I compute integer exponentiation modulo a constant in C

Question

I want to compute a^b mod c where a, b and c are integer values. Is there an efficient way to do this?

pow(a, b) % c does not seem to work.

Answer 1

For this task, the pow function defined in <math.h> is not the right tool for multiple reasons:

• computing large numbers with the pow() function may overflow the magnitude representable by the type double and the precision will not provide the low order digits that are required for the modulus operation.
• depending on its implementation in the C library, pow() might produce non-integral values for integral arguments, whose conversion to int could truncate to incorrect values.
• the % operation is not defined on the double type.
• to avoid losing the low order bits, one should perform the modulus operation at every step.
• for a test, it is not what the examiner expects.

One should instead compute the power and modulus iteratively in a loop:

unsigned exp_mod(unsigned a, unsigned b, unsigned c) {
    unsigned p = 1 % c;
    while (b-- > 0) {
        p = (unsigned long long)p * a % c;
    }
    return p;
}

For very large values of b, iterating will take a long time. Here is an efficient exponentiation algorithm that reduces the time complexity from O(b) to O(log b):

unsigned exp_mod_fast(unsigned a, unsigned b, unsigned c)) {
    unsigned p;

    for (p = 1 % c; b > 0; b = b / 2) {
        if (b % 2 != 0)
            p = (unsigned long long)p * a % c;
        a = (unsigned long long)a * a % c;
    }
    return p;
}

As suggested by rici, using type unsigned long long for the intermediary product avoids magnitude issues for large values of a. The above functions should produce the correct results for all 32-bit values of a, b and c, except for c == 0.

The initial step p = 1 % c is necessary to produce the result 0 for c == 1 && b == 0. An explicit test if (c <= 1) return 0; might be more readable and would avoid the undefined behavior on c == 0. Here is a final version with a test for special cases and one less step:

unsigned exp_mod_fast(unsigned a, unsigned b, unsigned c)) {
    unsigned p;

    if (c <= 1) {
        /* return 0 for c == 1, which is the correct result */
        /* also return 0 for c == 0, by convention */
        return 0;
    }
    for (p = 1; b > 1; b = b / 2) {
        if (b % 2 != 0) {
            p = (unsigned long long)p * a % c;
        }
        a = (unsigned long long)a * a % c;
    }
    if (b != 0) {
        p = (unsigned long long)p * a % c;
    }
    return p;
}

A more general analysis is available on the Wikipedia article titled Modular exponentiation.