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?
-
4The cpython source repo is at hg.python.org/cpython– WoobleMar 9, 2011 at 14:07
-
2...under Objects/longobject.c:long_pow() (as JimB had already commented).– smciJul 19, 2011 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.
-
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. Aug 31, 2015 at 12:53
-
Note that
long_pow
is now defined at another line in that file: github.com/python/cpython/blob/master/Objects/…– JohanCDec 8, 2019 at 21:57 -
@JohanC I've updated the link to include the commit hash, so it doesn't get out of date anymore. Dec 9, 2019 at 8: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;
}
-
@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) May 7, 2013 at 2:24
-
@stackuser: It appears to be working fine in the following demonstration: ideone.com/sYzqZN May 7, 2013 at 12:50
-
5Can anyone explain why this solution works? I am having trouble understanding the logic behind this algorithm. May 24, 2015 at 3: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
– FabianoJun 8, 2016 at 3: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. Jun 8, 2016 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.
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.
-
7
math.pow()
does't have the modulo argument, and isn't the same function as the builtinpow()
. Also FYI, Psyco is getting pretty stale, and no 64-bit support. NumPy is great for serious math.– JimBMar 9, 2011 at 14:30
Python uses C math libraries for general cases and its own logic for some of its concepts (such as infinity).
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