How did Python implement the built-in function pow()?

If a, b and c are integers, the implementation can be made more efficient by binary exponentiation and reducing modulo c in each step, including the first one (i.e. reducing a modulo c before you even start). This is what the implementation of long_pow() does indeed. The function has over two hundred lines of code, as it has to deal with reference counting, and it handles negative exponents and a whole bunch of special cases.

At its core, the idea of the algorithm is rather simple, though. Let's say we want to compute a ** b for positive integers a and b, and b has the binary digits b_i. Then we can write b as

b = b_0 + b1 * 2 + b2 * 2**2 + ... + b_k ** 2**k

ans a ** b as

a ** b = a**b0 * (a**2)**b1 * (a**2**2)**b2 * ... * (a**2**k)**b_k

Each factor in this product is of the form (a**2**i)**b_i. If b_i is zero, we can simply omit the factor. If b_i is 1, the factor is equal to a**2**i, and these powers can be computed for all i by repeatedly squaring a. Overall, we need to square and multiply k times, where k is the number of binary digits of b.

As mentioned above, for pow(a, b, c) we can reduce modulo c in each step, both after squaring and after multiplying.

You might consider the following two implementations for computing (x ** y) % z quickly.

In Python:

def pow_mod(x, y, z):
    "Calculate (x ** y) % z efficiently."
    number = 1
    while y:
        if y & 1:
            number = number * x % z
        y >>= 1
        x = x * x % z
    return number

In C:

#include <stdio.h>

unsigned long pow_mod(unsigned short x, unsigned long y, unsigned short z)
{
    unsigned long number = 1;
    while (y)
    {
        if (y & 1)
            number = number * x % z;
        y >>= 1;
        x = (unsigned long)x * x % z;
    }
    return number;
}

int main()
{
    printf("%d\n", pow_mod(63437, 3935969939, 20628));
    return 0;
}