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?

- 34,358
- 48
- 134
- 179
-
4The cpython source repo is at http://hg.python.org/cpython/ – Wooble Mar 09 '11 at 14:07
-
2...under _Objects/longobject.c:long_pow()_ (as JimB had already commented). – smci Jul 19 '11 at 20:45
6 Answers
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.

- 574,206
- 118
- 941
- 841
-
1
-
2@BenSandler: Because _a_ ≡ _a'_ (mod _c_) and _b_ ≡ _b'_ (mod _c_) imply _ab_ ≡ _a'b'_ (mod _c_), or in other words, it doesn't matter whether you first reduce _a_ and _b_ modulo _c_ and then multiply them, or multiply them first and then reduce modulo _c_. See [the Wikipedia article on modular arithmetic](https://en.wikipedia.org/wiki/Modular_arithmetic#Integers_modulo_n). – Sven Marnach Aug 31 '15 at 12:53
-
Note that `long_pow` is now defined at another line in that file: https://github.com/python/cpython/blob/master/Objects/longobject.c#L4246 – JohanC Dec 08 '19 at 21:57
-
@JohanC I've updated the link to include the commit hash, so it doesn't get out of date anymore. – Sven Marnach Dec 09 '19 at 08:38
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;
}

- 21,433
- 16
- 79
- 117
-
@Noctis, I tried running your Python implementation and got this:TypeError: ufunc 'bitwise_and' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe'' ---- As I'm learning Python right now, I thought you might have an idea about this error (a search suggests it might be a bug but I'm thinking that there's a quick workaround) – stackuser May 07 '13 at 02:24
-
@stackuser: It appears to be working fine in the following demonstration: http://ideone.com/sYzqZN – Noctis Skytower May 07 '13 at 12:50
-
5Can anyone explain why this solution works? I am having trouble understanding the logic behind this algorithm. – kilojoules May 24 '15 at 03:59
-
2@NoctisSkytower, what would be the benefit of this considering the native python `pow()` builtin function supports this as well and seems faster? `>>> st_pow = 'pow(65537L, 767587L, 14971787L) >>> st_pow_mod = 'pow_mod(65537L, 767587L, 14971787L)' >>> timeit.timeit(st_pow) 4.510787010192871 >>> timeit.timeit(st_pow_mod, def_pow_mod) 10.135776996612549` – Fabiano Jun 08 '16 at 03:25
-
7@Fabiano My function is not supposed to be used. It is simply an explanation of how Python works behinds the scenes without referring to its source in C. I was trying to answer **wong2's** question about how `pow` was implimented. – Noctis Skytower Jun 08 '16 at 13:18
I don't know about python, but if you need fast powers, you can use exponentiation by squaring:
http://en.wikipedia.org/wiki/Exponentiation_by_squaring
It's a simple recursive method that uses the commutative property of exponents.

- 8,601
- 7
- 58
- 90
Implement pow(x,n) in Python
def myPow(x, n):
p = 1
if n<0:
x = 1/x
n = abs(n)
# Exponentiation by Squaring
while n:
if n%2:
p*= x
x*=x
n//=2
return p
Implement pow(x,n,m) in Python
def myPow(x,n,m):
p = 1
if n<0:
x = 1/x
n = abs(n)
while n:
if n%2:
p*= x%m
x*=x%m
n//=2
return p
Checkout this link for explanation

- 505
- 7
- 12
Python uses C math libraries for general cases and its own logic for some of its concepts (such as infinity).

- 4,638
- 3
- 38
- 49

- 2,852
- 2
- 24
- 34
Line 1426 of this file shows the Python code that implements math.pow, but basically it boils down to it calling the standard C library which probably has a highly optimized version of that function.
Python can be quite slow for intensive number-crunching, but Psyco can give you a quite speed boost, it won't be as good as C code calling the standard library though.

- 10,682
- 3
- 35
- 38
-
7`math.pow()` does't have the modulo argument, and isn't the same function as the builtin `pow()`. Also FYI, Psyco is getting pretty stale, and no 64-bit support. NumPy is great for serious math. – JimB Mar 09 '11 at 14:30