I have to write a program to calculate a**b % c
where b
and c
are both very large numbers. If I just use a**b % c
, it's really slow. Then I found that the built-in function pow()
can do this really fast by calling pow(a, b, c)
.
I'm curious to know how does Python implement this? Or where could I find the source code file that implement this function?
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.